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Recent experimental achievements in controlling ultracold gases in optical lattices open a new perspective 
on quantum many-body physics. In these experimental setups it is possible to study coherent time evolution 
of isolated quantum systems. These dynamics reveal new physics beyond the low-energy properties usually 
relevant in solid-state many-body systems. In this paper we study the time evolution of antiferromagnetic 
order in the Heisenberg chain after a sudden change of the anisotropy parameter, using various numerical and 
analytical methods. As a generic result we find that the order parameter, which can show oscillatory or non- 
oscillatory dynamics, decays exponentially except for the effectively non-interacting case of the XX limit. For 
weakly ordered initial states we also find evidence for an algebraic correction to the exponential law. The study 
is based on numerical simulations using a numerical matrix product method for infinite system sizes (iMPS), 
for which we provide a detailed description and an error analysis. Additionally, we investigate in detail the 
exactly solvable XX limit. These results are compared to approximative analytical approaches including an 
effective description by the XZ-model as well as by mean-field, Luttinger-liquid and sine-Gordon theories. This 
reveals which aspects of non-equilibrium dynamics can as in equilibrium be described by low-energy theories 
and which are the novel phenomena specific to quantum quench dynamics. The relevance of the energetically 
high part of the spectrum is illustrated by means of a full numerical diagonalization of the Hamiltonian. 
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I. INTRODUCTION 

A. Motivation 

Developing a profound understanding of interacting quan- 
tum many-body systems is one of the important challenges 
in modern physics. Potential applications of correlated quan- 
tum systems have driven decades of theoretical and experi- 
mental investigation. There are very well understood classes 
of many-body systems which essentially behave as ensembles 
of non-interacting particles, with Landau's Fermi liquid the- 
ory as the most prominent example. This picture can how- 
ever break down in the presence of strong correlations or in 
reduced dimensions. For high-temperature superconductiv- 
ity [1], in the quantum-Hall effect or for transport in semi- 
conductor nanodevices interactions lead to intricate quantum 
many-body phenomena for which the theoretical basis is still 
incomplete. 

A variety of analytical and numerical techniques have been 
developed in order to study models of correlated systems. 
However, the degrees of freedom in quantum mechanics grow 
in general exponentially with the number of particles - a fun- 
damental problem which can make the analytical approach 
highly complex and restricts the applicability of computer 
simulations. An alternative to the analytical and numerical 
treatment of the quantum many-body problem has been pro- 
posed by R. Feynman [2], who introduced the idea of quan- 
tum simulation: instead of solving the highly complex theory 
on a computer, the model could be tested directly by means 
of an artificially engineered quantum system. Over the last 
decade, remarkable experimental setups have been developed 
to produce and control dilute ultracold atomic gases in optical 
lattices [3]. These systems are very promising candidates for 
the realization of Feynman's idea of quantum simulation. The 
key to reach collective quantum states of atomic gases was 
the development of laser and evaporative cooling techniques, 
which allow to go down to temperatures in the nanokelvin 
range and led to the first realizations of Bose-Einstein con- 
densation in 1995 [4-6]. Subsequently, using optical lattices, 
it became possible to drive ultracold gases into bosonic [7, 8] 
and fermionic [9-11] correlated states. 

A particularity of ultracold atomic systems is the controlla- 
bility of internal parameters, which relies on the development 
of magnetic and optical traps of various geometries [3], and 
the tuning of atom-atom interactions by means of Feshbach 
resonance in an external magnetic field [12, 13]. These pa- 
rameters can be changed quickly, producing a so-called quan- 
tum quench, which allows to generate various types of non- 
equilibrium situations [14—20]. Unlike solids, where the elec- 
tronic system suffers from dissipation due to the coupling to 
lattice phonons, atomic gases are almost perfectly isolated 
quantum systems and exhibit coherent dynamics over large 
periods of time. Since the coherent dynamical processes are 
determined exclusively by the intrinsic properties of system, 
it allows to investigate specifically the non-linear interaction 
effects. This is a unique situation, not available in usual solid- 
state experiments, where dynamical effects beyond linear re- 
sponse are in general intricate. The theoretical study of the 



non-perturbative many-body aspects of non-equilibrium dy- 
namics is the main topic of this paper. We namely focus on 
quantum spin chains, for which corresponding experiments 
are currently under development [19, 21]. 

Due to weak dissipation in ultracold atomic gases non- 
equilibrium properties are important even if equilibrium as- 
pects of these systems shall be investigated. When attempting 
to use ultracold gases as quantum simulators for certain equi- 
librium problem, one usually prepares the system in an un- 
corrected initial state, e.g. a Bose-Einstein condensate, and 
drives it into a correlated state by a slow change of parame- 
ter [7, 22]. However, it has been found that in the vicinity of a 
phase transition the time scales needed for a correlated state to 
equilibrate can become exceedingly large [23-28]. It is there- 
fore mandatory to establish relations between equilibrium and 
non-equilibrium systems to overcome this problem. 

The relevance of non-equilibrium dynamics of cold atoms 
goes beyond the investigation of fundamental questions of 
quantum statistical mechanics. There are possible practical 
applications in quantum metrology [29] and quantum infor- 
mation processing [30-33] and relations to open questions in 
cosmology [34, 35]. However, before the ultracold quantum 
gases can be routinely applied in the context of such problems, 
a number of open experimental challenges need to be solved. 
By providing exact results for realistic experimental setups in 
this article we intend to support ongoing efforts in improving 
the control of ultracold atomic gases. 

We will study the emerging dynamics of the order parame- 
ter of a XXZ Heisenberg chain prepared in the classical (un- 
correlated) Neel state, which can be realized in experiment, 
but, in order to get a deeper insight into the problem, general 
antiferromagnetic initial states are also considered. Our spe- 
cial interest concerns the effect of the quantum phase transi- 
tion crossed when tuning the magnetic anisotropy parameter. 

Exact results based on numerical calculations are pre- 
sented. Furthermore, alternative approximative approaches 
are applied. The applicability of the analytical tools, which 
have been very successful in the description of equilibrium 
phenomena, turns out to be strongly restricted for the non- 
equilibrium problem under consideration. We identify the ap- 
parent problems in the standard approximations and point out 
in which direction these approaches should be extended in or- 
der to capture the main features of the quantum quench dy- 
namics. 



B. Brief review on non-equilibrium dynamics 

Mostly in relation to transport phenomena, non-equilibrium 
problems have been subject to intensive theoretical investi- 
gation over many years (e.g. Ref. [36]). However, non- 
equilibrium transport, which can be seen as a result of per- 
turbations (voltage biases) at the edges of the system, is fun- 
damentally different from quench dynamics, where the pa- 
rameter change is global. More closely related to a quantum 
quench are highly excited electronic states in solids, gener- 
ated in femtosecond pump-probe spectroscopy [37-39]. Nev- 
ertheless, in these systems decoherence times are short and the 
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dissipative processes strongly contribute to the emerging dy- 
namics. Consequently, concepts developed for transport phe- 
nomena and dynamics in condensed matter systems are not 
necessarily appropriate to quenches in ultracold atomic sys- 
tems. Except for pioneering works on quench dynamics in 
the 1970's [40-43], specific theoretical research on quench 
dynamics has only started in recent years, stimulated by the 
experimental developments in ultracold atomic physics. In 
these works, which shall be briefly summarized in this sec- 
tion, two main lines have been followed. A first aspect is the 
study of the nature of the quasi-stationary states in the long- 
time limit. As demonstrated in an experiment of Kinoshita 
et al. [16], these non-equilibrium states can exhibit striking 
properties for specific types of interactions. Another approach 
explicitly focuses on the characteristics of the time evolution 
after the quench - experimental examples are the oscillations 
[14] or the dephasing [17] of the superfluid phase. It turns out 
to be an ambitious challenge to establish relations between dy- 
namical phenomena and the details of the microscopic model, 
such as integrability and dimensionality. Although numerous 
remarkable theoretical efforts revealed a number of interesting 
phenomena, many aspects of relaxation dynamics and equili- 
bration, which shall be discussed in detail in this work, remain 
unclear. 

The effective description of many-body systems by means 
of low-energy theories, captured within the renormalization 
group framework [44], has proven sufficient for the theoret- 
ical understanding of a broad range of equilibrium phenom- 
ena. Therefore the application of renormalization group ideas 
to non-equilibrium dynamics seems to be a promising ap- 
proach. Along this way diagrammatic techniques [45^48] and 
the solutions of the dynamics of field-theoretical models at the 
renormalization group fixed point [49-57] were developed. 
Providing a generic view on the quench problem for critical 
theories, the work of Calabrese and Cardy [52, 53] based on 
conformal field theory has to be emphasized. While for con- 
tinuum systems field-theoretical models were successfully ap- 
plied to generic quantum quenches [17, 54], it has to be clari- 
fied under what conditions they all provide accurate decription 
of lattice systems. Unclear is also the range of applicability of 
semiclassical theories [58-60]. 

For a restricted class of problems the time-evolution can be 
calculated exactly, e.g. for Jordan-Wigner diagonalizable XY- 
chains [40-43, 61-64] or the i-Hubbard-chain [65]. A major 
drawback of these exactly solvable models is that the possibil- 
ity of their representation in terms of non-interacting particles 
apparently leads to very specific relaxation phenomena, which 
are not generic not only for non-integrable, but also for more 
complicated integrable models. For instance, it is question- 
able whether the generalized Gibbs ensemble, which has been 
proposed for the description of quasi-stationary states of in- 
tegrable models [63], is a useful concept beyond the simple 
Jordan-Wigner diagonalizable cases [64, 66]. For the more 
general Bethe-ansatz solvable models, it has not yet been pos- 
sible to extract dynamics, except for the Richardson [67] and 
the Lieb-Liniger models [68]. 

In view of the high complexity of the quench dynamics, 
efficient unbiased numerical approaches are crucial to gain 



deeper insight. Using exact diagonalization [69-71] it is pos- 
sible to calculate the dynamics of small systems over exceed- 
ingly long times. For larger (but one-dimensional) systems the 
density matrix renormalization group (DMRG) [72-75] can 
be applied. Although only for finite times, dynamics of spin- 
chains (respectively spinless fermions) [76-82] and bosonic 
lattice models [70, 83-85] have been evaluated. Recently, the 
dynamical mean field theory has been applied to fermionic 
models in the limit of infinite dimensions [86-89]. 

C. Basic setup and general discussion 

The Heisenberg model is a paradigm in the theory of mag- 
netism and strongly correlated systems in general. In ap- 
pendix A we derive how the model can be realized with ultra- 
cold two-level atoms in various geometries of optical lattices. 
For instance, it is possible generate a one-dimensional XXZ 
Heisenberg model, 

H = J E i s PU + s * s hi + A ^ + i) - "> 
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where the sign and the strength of the exchange coupling J 
and A can be tuned dynamically. The XXZ model is in- 
tegrable and its eigenstates can be constructed by the Bethe 
ansatz. In the case of antiferromagnetic couplings J > 0, the 
anisotropy parameter triggers a quantum phase transition from 
a gapless "Luttinger liquid" phase (0 < A < 1) to a gapped, 
Ising-ordered antiferromagnetic phase (A > 1). The main 
features of the model at equilibrium and its field-theoretical 
formulation are given in appendix B. 

The non-equilibrium dynamics (1) shall be investigated in 
the following quantum quench: At time t < the system 
is prepared in a ground state \i/jq) with long-range antiferro- 
magnetic order. The corresponding anisotropy parameter is 
denoted as Ao, Ao > 1. Among the aniferromagnetic equi- 
librium states the Neel state, 

|V)N&i = |UT---ITi), (2) 

which corresponds to the limit Ao^oo, has already been re- 
alized in experiment [19] and will attract our special attention. 
At t = the system is pushed out of equilibrium by changing 
the strength of the interaction, A < Ao, and the dynamics 
emerging at t > are studied. 

In the context of optical lattices, where the system is well 
isolated and no phonons are present, dissipation can be ne- 
glected in a first approximation. Also, being interested in 
quantum effects, we set T = 0. Finite temperature may be- 
come relevant for weak magnetic exchange interactions in the 
ultracold atomic setup, but how to investigate efficiently the 
non-equilibrium problem at T > is still an unsolved prob- 
lem. Under these assumptions, the dynamics is formally de- 
scribed by the solution of the Schrodinger equation, 

Mi)) = e- mt \4>o) • (3) 

we set h = 1 throughout this paper. Involving a priori all 
the energy scales of the many-body Hamiltonian, the calcula- 
tion of the time evolution of the wave function (3) is highly 
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complex. When approaching the problem analytically, one 
is forced to intoduce an appropriate approximation - the ad- 
vantages and drawbacks of various approaches will be inves- 
tigated in this work. When using numerics the dynamics (3) 
can be solved by fully diagonalizing the Hamiltonian H. In 
section VIII we apply the full numerical diagonalization ap- 
proach. Highly efficient routines have been developed for this 
purpose [90], which can nevertheless be used only for small 
system sizes (up to 20 lattice sites). More efficient and ap- 
plicable directly in the thermodynamic limit are matrix prod- 
uct states (MPS), which will be used for the simulation of the 
general quench dynamics of the XXZ model. For a detailed 
description of the MPS method see appendix C. 

To describe the dynamics of the state \tp(t)) we mainly fo- 
cus on the antiferromagnetic order parameter, 

ms(t) = ^(-iy{S^(t)). (4) 

3 

Since the state \ipo) is invariant under translation and subse- 
quent spin-inversion, ±m s (t) corresponds to the local mag- 
netization at any site of the lattice. It will also be useful to 
look at the frequency distribution f ms (e), which resolves the 
contributions to the dynamics in energy space, 

m s (t) = J dee~^f ms (e). (5) 

The staggered magnetization is not only the natural observ- 
able characterizing the ordering of antiferromagnetic states, 
but also reflects the properties of the local density matrix of a 
single site. For describing non-local properties we choose the 
equal-time connected spin-spin correlation function, 

G?(i,t) = ^£<^(t)S?(*)> ~ (S? +t (t))(S?(t)) . (6) 

z 

Before going into the study of the many-body dynamics of 
the Hamiltonian (1), it is worthwhile considering the case of 
only 2 spins. A corresponding experiment has been carried out 
by Trotzky et al. [19] by loading 87 Rb atoms in the hyperfine 
states | i) = \F = l,m F = -1) , | t) = \F = 1, m F = 1} , 
into an array of double-well potentials. The initial Neel state 
was generated using a magnetic field gradient transferring the 
effective spins in each double well from a triplet-bond state 
into an antiferromagnetic one, IV'o) = |TD- The dynamics are 
in this special case independent of A and can be described as 
Rabi oscillations between ||J.) and |J.f) states, 

\m) = cos(Jt/2)|Ti) + * sin( Jt/2)||T) ■ (7) 

Hence, the antiferromagnetic order parameter descibes an os- 
cillatory behaviour, m s (t) = \ cos(Jt), where the Rabi fre- 
quency is set by the exchange coupling J, which was indeed 
observed in the experiment [19]. 

Although, as we shall see in section II, in a many-body 
system such Rabi-like oscillations may survive, the dynam- 
ics become much more intricate when going to large system 
sizes. On a heuristic level the initial state may be regarded as 



a bunch of excitations of the Hamiltonian H, whose dynam- 
ics gives rise to the propagation of correlations throughout the 
system. For spin models with sufficiently local interactions, 
Lieb and Robinson [91] have proven that this propagation 
takes place within a light-cone - the deviation of a correlation 
function from its initial value becomes exponentially small for 
distances £ > 2ut, where u is the maximum velocity of exci- 
tations in the system. For an isolated but arbitrarily large sys- 
tem this means that relaxation to a stationary state can only 
be observed for subsystems of size £ < 2ut. This light-cone 
effect has been more precisely described in the framework of 
boundary conformal field theory [92], which predicts an expo- 
nential decay of the correlations in the long-time limit. These 
short-range correlations are in contrast with the entanglement 
properties of the non-equilibrium problem. It has been shown 
[92] that the entanglement entropy of a subsystem of size £ 
grows linearly with time if 2ut < £ and saturates to a value 
proportional to £ if 2ut > £. 

It is an open question, under what conditions the station- 
ary state in the long-time limit can be described by a statis- 
tical ensemble at a finite temperature, meaning that thermal- 
ization occurs. There are several examples for which this is 
not the case [63, 71, 78, 79, 83], and the extended Gibbs 
ensemble [63], which takes into account the constraints of 
the non-dissipative dynamics, or the micro-canonical ensem- 
ble [66, 69, 93] are possible candidates for describing steady 
states. Whether the integrability is a necessary condition for 
the absence of thermalization remains unclear. It has been 
pointed out that the absence of thermalization can be associ- 
ated with a non-perturbative behavior, which is not related to 
the integrability of the underlying Hamiltonian [71]. 

Here, it will be shown that in the long-time limit the an- 
tiferromagnetic order vanishes in all cases, hence, at least 
for this local quantity, thermalization is observed - in a one- 
dimensional system no long-range order is possible at finite 
temperatures. This does not necessarily imply thermalization 
for correlation functions. Indeed, in section III we present re- 
sults which indicate the absence of thermalization in the spin- 
spin correlations (6). However, the correlation functions ex- 
hibit somewhat slow relaxation dynamics and it is difficult to 
extract steady-state properties from the rather short accessible 
times that can be achieved numerically. 

Nevertheless, interesting dynamical effects are present also 
at short times. Their characterization as a function of the ini- 
tial state and the interaction parameter will be investigated. 
The magnetic order parameter turns out to be a good observ- 
able for the quantitative extraction of non-trivial time scales. 
Here, where the initial state can be characterized by the gap 
parameter A s (more precisely, the inverse correlation length), 
one expects that the typical time scale of the relaxation dy- 
namics is given by A J 1 and the length scales, which depend 
on the momentum distributions in the initial states, should be 
of the order of u/A s , where u is given by the velocity of 
quasi-particles (spin-waves). 

In the solution of the quench dynamics for conformally in- 
variant theories [94] of Calabrese and Cardy [52, 53] these 
qualitative arguments were put on a solid ground: The ini- 
tial state enters into the framework of quantum field theory as 
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TABLE I: Exact analytical and numerical results for the quench in the XXZ model. See sections II - VI for details. 



initial coupling 
state 


asymptotic law 


relaxation times 
(if applicable) 


Sec. II: Exact analytical calculation in the XX limit 


Neel A = 

SDW A = W 


t~i cos(2 Jt - f ) 
^-i e - 2Ast + fyt-i cos(2Jt - f ) 


T\ « , r 2 — > oo 

T l = 2^7 ' T 2 00 


Sec. Ill: Numerical iMPS calculation of the XXZ model 



Neel A > 
Neel A > 1 
A > 1 A = 
Ao > 1 A > 
A > 1 A > A > 1 


e~ t/T2 cos(u;t + (j>) b 
e -t/n 

C\t~h- t/T1 +C 2 t~? cos(wf + </>) * 
dt-h~ t/T1 +C 2 e- t/T2 cos(Lut + (l>) b 
e -t/n 


ri » , r 2 ~ log A 
n ~ A 2 

n ~ ,t 2 ~ log A 

- 1 1 1 |- 2 


Sec. IV: Mean field theory 


Neel 1 > A > 
Neel A > 1 


t~3 [d cos(wii + 0i) + C 2 cos(u! 2 i + <fe)} 
const. 


n « , r 2 — > 00 


Sec. V: XZ model 


Neel A > 1 
Neel A < 1 


e -*/n 

e~ i/T2 (cos 2 (o;t) - const) 


Tl ~ A 2 

T 2 ~ A -1 


Sec. VI: Luttinger model 


KG A > 


e -*/ri 


2 

Tl — A'jrA,, 



"Valid in an intermediate time regime (See section III). 

6 Only rough correspondence, deviations of the order of 30% are possible. 



a finite slab width, r e , the extrapolation length which stands 
for the renormalization-group distance of the initial state from 
the fixed point of the gapped theory [95]. To first order, this is 
given by the inverse gap, here r e ~ A J 1 . Using a conformal 
transformation, the slab geometry is mapped onto the semi- 
infinite plane, for which, by means of boundary conformal 
field theory [96], the properties of the correlation functions 
can be extracted. 

The results of Calabrese and Cardy [53] do apply to the 
quench in the XXZ model if the discussion is restricted to the 
low-energy modes in the gapless regime |A| < 1, here cap- 
tured by the Luttinger model [see appendix B, Eq. (B9)]. For 
the staggered magnetization as a local observable the outcome 
is 

m s {t) ~ e'^i , (8) 

where r e ~ A J/ 1 . 

However, several remarks concerning the applicability of 
the conformal field theory results to the quench in the XXZ 
chain are in place. First, the initial state is treated on a pertur- 
bative level in terms of a renormalization-group distance from 
the fixed point and simply characterized by the gap parame- 
ter. It is questionable whether in this framework it is possible 
to correctly take into account the physics of the antiferromag- 
netic states, especially those close to the critical point (i.e. far 
from the antiferromagnetic fixed point), described by the sine- 
Gordon model. Second, within the field theory it is impossi- 
ble to treat lattice effects, which are expected to emerge if the 



energy of the quasi-particles forming the initial state is of the 
order of the bandwidth A - a situation which is realized for in- 
stance by the Neel state (2). As a simple example of a lattice 
effect we presented the Rabi-oscillations in the two-spin sys- 
tem (7), with the frequency set by the magnetic exchange J. 
Macroscopic order parameter oscillations following a quan- 
tum quench have been predicted to appear in a variety of sys- 
tems [58, 62, 97-99]. In this work we will characterize Rabi- 
like oscillations and investigate origins of dephasing in the 
presence of many-body correlations. A particular property of 
the quench in the XXZ chain illustrates the novel aspect of the 
non-equilibrium dynamics in many-body lattice models: The 
time-evolution of m s (t) is invariant under the change of sign 
A — ► — A. Ferro- and antiferromagnetic Hamiltonians exhibit 
identical dynamics despite their completely different elemen- 
tary excitations. As a third point restricting the applicability 
of the conformal field theory result, we mention that a confor- 
mal theory does not capture the case of a parameter quench 
into the gapped phase, A > 1. Here this regime will be ad- 
dressed using a sine-Gordon description of the XXZ model. 



D. Summary of the results 

The further content of the paper is organized as follows: 
The non-equilibrium dynamics in the XX limit of the Heisen- 
berg model, which can be solved in a simple way by means of 
Jordan- Wigner transformation, is analyzed in section II. Nu- 
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merical results for the general case are given in section III and 
approximative approaches in sections IV-VI. In section VIII 
an exact diagonalization analysis of the spectrum of the XXZ 
model is carried out before presenting the conclusions. In ap- 
pendix A we describe the experimental realization of quan- 
tum magnetic systems in optical lattices. The well-established 
properties of antiferromagnetic states and equilibrium phase 
transitions in the context of the Heisenberg model in one di- 
mension are reviewed in appendix B. The description and an 
error analysis of the matrix product algorithm is provided in 
appendix C. 

Our results for the non-equilibrium dynamics of the stag- 
gered magnetization are summarized in Table I. We find es- 
sentially two types of relaxation dynamics: non-oscillatory 
dynamics, characterized by a relaxation time n, and oscilla- 
tory dynamics with a frequency u> and an associated relaxation 
time T2- An important result is that for non-zero A we find a 
fundamentally new mode of many-body dynamics which al- 
ways leads to exponential decay of the staggered moment re- 
gardless of whether the short-time dynamics is oscillatory or 
not. In contrast with the oscillation frequency, which is set by 
the exchange interaction, the relaxation time is an emergent 
scale generated by the highly correlated dynamics and hence 
cannot be simply related to the microscopic parameters. We 
find divergent relaxation times, ti^oo in the limit A — > 
and T2^oo if A — > oo. For the particular case of the Neel 
state, we find that the relaxation times essentially vanish in 
the vicinity of the critical point, A > 1. 

Table I also shows to what extent approximative methods, 
which take into consideration only a particular aspect of the 
Hamiltonian, are applicable to the non-equilibrium problem. 
The mean-field approximation for example leads to contradic- 
tions with the unbiased numerical results - an algebraic decay 
for A < 1 and a non-vanishing asymptotic value of the stag- 
gered moment for A > 1 [97]. In the case of the initial Neel 
state, comparing the low-energy result of conformal field the- 
ory with the numerics, the immediate relaxation n ~ is in 
principle in agreement with A s — >oo in Eq. (8). However, the 
oscillations dominate the long-time dynamics, and are, as ex- 
pounded before, not captured by the field theory. If the initial 
state is close to the critical point, an exponential relaxation 
similar to Eq. (8) is found, however, an additional algebraic 
prefactor appears to be present. In our treatment of the Lut- 
tinger model this effect is also not seen, but the results from 
conformal field theory (8) are reproduced. 



A. Initial Neel state, Ao = oo 

In the fermionic picture, the Neel state reads as a charge 
density wave, 



NAo>= n (4 



-k+ 



J|o). 



(9) 



: -<k<i 



The fermionic operators are easily represented in the Heisen- 
berg picture, 



Jtekclck 



c k e 



-ite k cick 



c k e 



(10) 



Hence, the dynamics of the XX chain prepared in the Neel 
state, in analogy with the two-site model (7), takes the form of 
Rabi oscillations between charge-density waves with different 
sublattice magnetizations, 



W)) 



n (< 



MktJ 



c k + 6 



4+Jlo) 



(ii) 



-f <k<l 



The relaxation of the staggered magnetization can be seen as 
a dephasing process, driven by inhomogeneous Rabi frequen- 
cies in fc-space, 



.(*) = ^ E e^iMcUk+M 

1 E -< 



J2t k t 



N _ 

2 <k< 2 

In the thermodynamic limit, 

1 1 

= - / dkcos(2te k ) = -J (2Jt) , 
* Jo 2 



(12) 



(13) 



where Jo denotes the zeroth Bessel function of the first kind. 
The underlying frequency distribution (5) ranges over a band 
of width 4 J, 



/ m ,(e) =6(2J-e)9(e + 2J) 



1 



V4 J 2 - e 2 



(14) 



8(e) being the Heaviside function. High-energy modes with a 
vanishing velocity at the band edge, \e k \ = J, dominate the 
long-time limit of (13) and give rise to the oscillations with a 
frequency set by the bandwidth, 



H. XX MODEL, A = 

It is particularly illustrative to study the exactly solvable 
case of zero anisotropy (A = 0), where the Heisenberg 
Hamiltonian (1) can be represented in terms of free spinless 
fermions with a cosine dispersion relation (B4). For free 
fermions the non-equilibrium dynamics can be solved analyt- 
ically [100]. We study two cases: first, the Neel state as the 
initial condition, second, the case of the initial spin-density- 
wave state. 



m s (t) 



Jt»i 



1 / . 

cos(2 Jt ) 

4tt Jt y 4 ; 



(15) 



The exponent of the decay is a consequence of the 
quadratic dispersion at k = 0. In the XX limit it is also possi- 
ble to express the correlation function, G z c z (£,t), in terms of 
Bessel functions, 



G zz (£,t) 



Se,o - 1 



-in 



Tl/ 2 



-it/2 



dk cos(k£) cos(2te k ) 



(S e ,o- Jt(2Jt)) 



(16) 
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This results in slowly decaying, spatially oscillating correla- 
tions, 



vicinity of the band edges, 



G' c '{e,t) 



2nJt 



cos 2 (2Ji- £7r/2-7r/4) (17) 



Fig. 1 shows how the correlations evolve within the light cone 
£ < 2t. The magnitude of the wave-front decays as a power 
law in time. The negative sign reflects spinon characteristics 
[101] of the propagating correlations. 
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FIG. 1 : Correlation functions in the XX limit. Comparison of results 
obtained for full (a) and linearized (b) spectra. For the linearized 
spectrum we set A = 2 J as the effective bandwidth. 

Although it is possible to carry out the analysis of the XX 
model without any approximation, it is useful to investigate 
the result of restriction to a particular part of the spectrum. 
This provides information on the range of applicability of low- 
energy theories, which are candidates for treating the more 
complicated case of interacting systems. 

In the case of the linearized theory [appendix B, Eq. (B8)], 
the dynamics of the magnetization is characterized by oscilla- 
tions with a 1/t decay and cutoff-dependent period, 



m s {t) 



1 

At 



sin(At) 



(18) 



The cutoff A gives the correct periodic behavior if it is equal 
to the bare bandwidth (A = 2 J). The oscillatory behavior, a 
consequence of the presence of the lattice, is indeed not cap- 
tured in the continuum limit A/ J — > oo, where the oscilla- 
tions disappear. The power-law decay appears in the linear 
approximation being independent of the cutoff, but the expo- 
nent is overestimated by a factor of two compared to the case 
of the full dispersion. The energy distribution corresponding 
to the magnetization (18) is simply flat, 



f m ,(e) = 0(2A-e)0(2A + e) 



(19) 



A seemingly (in the context of equilibrium theories) un- 
conventional approach is the development of the modes in the 



Hq= J2 - J ( 1 - k2 ){ C l C k-4+n C k+^} 



(20) 



Jk 2 <A 



In the present case of non-equilibrium dynamics, we find, 
however, that the corresponding energy distribution, 



/m s 0) - 6{e + 2J)6{K-2J - e) 
+ 9(2 J - e)0(A - 2 J + e) 



1 



V2J+1 
1 

V2J-e ' 



(21) 
(22) 



provides the correct long-time limit if the cutoff is sufficiently 
large, 



1 7T 

Jt»A-! y/t 4 



(23) 



We now clearly understand the mechanism behind the de- 
phasing process in the free-fermion models: Rabi oscillations 
are present if there is a sharp step at the edge of the band. The 
dephasing of the oscillations is algebraic, t~ a , a = 1 if the 
frequency distribution is homogeneous and a = | in the case 
of the quadratic dispersion at the band edge. For the long- 
time behavior it is sufficient to stick to the modes at the edge 
of the band, the low-frequency part is effective only at short 
times t ~ J -1 . The reason for such behavior is best illus- 
trated in the analysis of the correlation functions for the linear 
spectrum. The result, as shown in Fig. 1, is a single coherent 
spinon mode traveling the light cone \2t — £\ = 0. For the 
staggered magnetization as a local observable this means that 
it relaxes as soon as the spinon mode moves over more than 
one lattice distance 2t > 1. In contrast to the case of the full 
dispersion, there are no oscillations within the light cone. We 
note that this immediate decay is in agreement with the result 
of conformal field theory (8), which predicts zero relaxation 
time for the Neel state due to its vanishing correlation length 
(inverse gap). 



B. Initial spin-density wave 

As an introduction to our discussion of quenches from cor- 
related antiferromagnetic states (i.e. quenches with 1 < Ao < 
oo), we consider the time evolution of weakly antiferromag- 
netic spin-density-wave states under the XX Hamiltonian [See 
appendix B, Eqs. (B4) and (B20)]. This section will provide 
a benchmark for the numerical results in section III and also 
discusses the applicability of effective low-energy theories to 
this quench. 

The time evolution of the staggered magnetization m s (t) 
in the XX model starting from a SDW state at t = (\ipo) = 
Il^/2<k<7r/2( u kcl + Vkc\ + J\0)) is determined by the co- 
efficients Uk and Vk, 

1 " 



fc=-i 
*■ dk 



— e UkVk , 
Ztt 



(24) 
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where we have taken the thermodynamic limit in the last equa- 
tion. With the coefficients obeying u^Vk 



dephasing process in energy representation reads 



the 



1 f° 

i s (t) = - \ 



cos(2ie) 



de 



2^ 



(25) 



For a weak SDW state (A s <C 1) there are two main con- 
tributions to the integral in Eq. (25). The first comes from 
the Fermi points e = 0, whereas the second originates in the 
square root singularities at e = ± J. We write these two con- 
tributions separately, 



m s {t) ^^-K (2A s t) + J (2Jt) 
it J 2 J 



(26) 




-2A.t 



J 



cos(2Ji-7r/4) 



In comparision with the case of the initial Neel state, in ad- 
dition to identical algebraically decaying oscillations (13) a 
non-oscillatory decay stemming from the low-energy part of 
the spectrum is obtained. This exponential behavior with 
an algebraic prefactor is characterized by the relaxation time 
t = (2A S ) _1 . Hence, for t > Aj 1 ln(J/A s ) the oscilla- 
tions on top of the non-oscillatory decay dominate the order- 
parameter dynamics. Nevertheless, unlike the case of the ini- 
tial Neel state, the low-energy modes contribute to the non- 
equilibrium dynamics over significant periods of time. 



III. INTERACTION QUENCH IN THE XXZ-MODEL - 
NUMERICAL STUDY 

In this section we first study the quench in the XXZ model 
starting from the Neel state. Subsequently ground states of the 
XXZ models at finite A = Ao will be considered. 

Unlike for A = 0, the problem is no longer analytically 
treatable and we have to resort to numerical techniques. In 
the iMPS algorithm (appendix C) we use 2000 states and a 
second-order Suzuki-Trotter decomposition with a time step 
5 — 10~ 3 J -1 for large A and up to 7000 states with S — 
10~ 2 J -1 for small A. An intermediate time regime Jt < 16 
can be reached, which exceeds in general greatly the short 
transient time. 



A. Initial Neel state, Ao = do 

An overview of the results for the initial Neel state is pre- 
sented in Fig. 2. For small anisotropies we find oscillations 
of the order parameter similar to those in the XX limit, but 
with the decay time decreasing upon approaching the isotropic 
point A = 1. In the easy-axis regime A > 1 of the XXZ 
model, the relaxation slows down again for increasing A and 
we observe non-oscillatory behavior for A ^> 1. 

Fig. 3 focuses on easy-plane anisotropy < A < 1. The 
results for < A < 0.4 are well described, for accessible 




FIG. 2: Dynamics of the staggered magnetization m a (t) in the XXZ 
chain initialized in a Neel state. Symbols correspond to numerical 
results, lines represent analytical results or fits by corresponding laws 
(27) and (28) (see text). 




0.001 



FIG. 3: Analysis of the decay of the oscillations in the XXZ model 
by plotting the absolute value of the staggered magnetization on a 
logarithmic scale. Symbols represent numerical results, solid curves 
correspond to fits by the function (27), straight lines point out the 
exponential decay. For anisotropies close to zero (A = 0.2, 0.4) an 
exponential law is obeyed for large periods of time in the numerically 
accessible time window, while for A = 0.6 the simulation shows an 
acceleration of the decay after a few oscillations. 



time scales, by exponentially decaying oscillations 

m s (t) oc e~'/ T2 cos(wi + cj>) . 



(27) 



The oscillation frequency is almost independent of the 
anisotropy, while the relaxation time t-i increases with de- 
creasing A. Logarithmic divergence of the relaxation time in 
the limit A — > is suggested by the fit shown in Fig. 4. The 
picture is less clear closer to the isotropic point. For the range 
0.5 < A < 1 there appears to be an additional time scale 
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FIG. 4: Relaxation times and oscillation period T = ^j- as a function 
of anisotropy in the XXZ model for the Neel initial state. Logarith- 
mic or algebraic laws are emphasized by solid lines. In the region 
close to the critical point of the XXZ model (indicated by the ques- 
tion mark) it becomes impossible to extract a relaxation time from 
the numerical results. 




FIG. 5: (a) Focus on the XXZ chain close to the critical point A = 1, 
where a steep decay of the initial magnetization is followed by a 
rather slow relaxation on tiny magnitudes which does not fit either 
of the generic behaviors (27,28). (b) Comparison of the XXZ chain 
(symbols) and the XZ chain (dashed lines) for strong anisotropies, 
solid lines correspond to an exponential fit. The dynamics of the stag- 
gered magnetization of the XXZ and XZ chains converge towards 
each other in the large-A limit. 



after which the oscillations start to decay even faster than ex- 
ponentially, simultaneously the period of the oscillations is 
reduced. Therefore, the relaxation times plotted in Fig. 4 are 
only valid within an intermediate time window, whose width 
shrinks upon approaching the critical point. 

For intermediate easy-axis anisotropies 1 < A < 3, the 
magnetization does not reach a stable regime within the nu- 




t-t Finite T, T=0.6 1 
t=4 



FIG. 6: The correlation functions obtained using iMPS for the initial 
Neel state. Symbols T denote quantum Monte Carlo results for the 
XXZ model at equilibrium at a temperature fixed by the energy of 
the non-equilibrium system. 



merically accessible time window [Fig. 5(a)]. The compli- 
cated behavior of m s (t) in this parameter range can be as- 
cribed to the interplay of processes at all energy scales. Never- 
theless, the numerical data suggest that the relaxation is fastest 
close to the isotropic point, in the range between A = 1 and 
A = 1.6. A simple generic type of behavior is recovered for 
large anisotropies A > 3. The numerical data in Fig. 5(b) in- 
dicates exponential relaxation of the staggered magnetization 



m s (t) oc e~ t/Tl . 



(28) 



The relaxation time scales roughly quadratically with A (Fig. 
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4). Oscillations do persist on top of the exponential decay, but 
they fade out quickly. 

We briefly describe the relaxation of the spin-spin correla- 
tion functions (6) as presented in Fig. 6. A more detailed 
study of these has been carried out by Manmana et al. [81]. 
For weak interactions (e.g. A = 0.6) the dynamics of corre- 
lation functions is still dominated by the spinon mode moving 
according to the light-cone [91, 92] set by the spin-wave ve- 
locity u [See appendix B, Eq. (Bl 1)], as it is the case at A = 
(Eq. (17), Fig. 1). For larger A, this mode is smeared off, in- 
stead, antiferromagnetic correlations build up. The strength of 
the short-range antiferromagnetic correlations increases as the 
anisotropy A is augmented. With the numerical method, how- 
ever, we are unable to reach sufficiently long times to calculate 
the quasi-stationary correlation length. It becomes neverthe- 
less clear that the correlations cannot be described in terms of 
a thermal ensemble. We evaluated the equilibrium correlation 
functions at a temperature corresponding to the energy of the 
system by means of quantum Monte Carlo simulations [176]. 
The resulting correlation functions depicted in Fig. 6 decay 
considerably faster than the non-equilibrium ones. 



B. Initial antiferromagnet, 1 < Ao < oo 

The Neel state is an entirely classical state with no quan- 
tum correlations. In order to generalize our results, we first 
study the case of small but finite correlations starting from the 
ground state for Ao = 4.0. We find that the picture gained 
from the initial Neel state remains qualitatively valid - the dy- 
namics of m s (t) is very similar to that in the case of the initial 
Neel state (Fig. 2). The corresponding relaxation times and 
periods are plotted in Fig. 7. For A close to zero, the behavior 
of t-2 is again close to a logarithmic law and the divergence of 
the relaxation time for A — > Ao goes like ti oc | — ^ | ~ 2 . 

We expect a qualitatively different behavior for a weakly 
ordered (more strongly correlated) initial state. In section II 
we have seen that for an initial spin-density-wave state and 
the XX Hamiltonian, in addition to the algebraically decay- 




FIG. 7: Relaxation times and oscillation period T = — as a func- 
tion of anisotropy in the XXZ model for the system prepared in the 
ground state for Ao = 4. Logarithmic or algebraic laws are empha- 
sized by solid lines. 



ing oscillations, an exponential relaxation exists, whose re- 
laxation rate is proportional to the gap of the initial state. In 
Fig. 8, where we show the results for the quench from an 
initial state with Ao = 1.5, oscillations are found on top of 
non-oscillatory relaxation. At A = 0, for sufficiently large t, 
the dynamics is similar to the SDW result (26), 



m s (t) ~ C 1 t~?e~ t/Tl + C 2 t~^ cos(ujt + 



(29) 



to very high accuracy, despite the fact that the spin-density 
wave is a different wave function than the ground state of the 
XXZ chain. The relaxation time t\ « 5.1J is slightly smaller 
than the one predicted by the SDW calculations, (2A S ) _1 k, 
5.8 J. The difference may be explained by the importance of 
short-range effects which are supposed to contribute to the 
non-equilibrium dynamics. As illustrated in Fig. 17, the cor- 
relations decay much faster at shorter distances than in the 
large distance asymptotics. 

For A > 0, in correspondence with the result for the ini- 
tial Neel state, we find that the oscillations are exponentially 
damped, while the non-oscillatory part remains qualitatively 
the same as in the XX limit, 



m s {t) ~ dt-h- 1 / 7 ' 1 + C 2 e~ t/T2 cos(ujt + , 



(30) 



In Fig. 9 we plot the fitting parameters for small A (0 < A < 
0.6) where formula (30) is well obeyed. t\ behaves similarly 
to (J/KA S ) - a law which is the natural extension of the 
non-interacting SDW result (26) to finite anisotropies using 
the same scaling as derived for Luttinger model [see section 
VI, Eq. (49)]. The logarithmic behavior of t 2 , apparent for 
A ^> 1, is not observed here. Oscillatory and non-oscillatory 
terms are superimposed. In the non-oscillatory term of (30) 
absence of the algebraic prefactor, as suggested by the field 
theoretical-result (8), can be clearly excluded on the basis of 
the numerical results. Pure exponential law (28) is however 
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FIG. 8: Dynamics of the staggered magnetization m s {t) in the XXZ 
chain prepared in an antiferromagnetic ground-state of the XXZ 
Hamiltonian with A = Ao = 1.5. Symbols correspond to numeri- 
cal results, lines represent analytical results or fits by corresponding 
laws (29), (30) and (28) (see text). 
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FIG. 9: Relaxation times and oscillation period T = — as a func- 
tion of anisotropy in the XXZ model for the system prepared in the 
ground state at Ao = 1.5. Solid lines are guides to the eye. ri is 
comparable to (J/KA S ). 
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found for A > 1. The intermediate regime 0.6 < A < 1 can 
not be described by either of the laws (28), (30). 



FIG. 10: Numerical solution of the mean field equations for the 
Hamiltonian (32) with the initial condition m s (0) = i. 



IV. MEAN FIELD 

Time-dependent mean field theory is one possibility to treat 
the dynamics of the XXZ model approximately. The mean- 
field approximation of the Hamiltonian (B2) at an instant of 
time t, is defined by expanding the interaction term to lin- 
ear order in fluctuations drij around the mean density, rij = 
(rij) + Srij, and by setting (rij) = 1/2 + (— iym s , 

7T 

H M F(t) = -J ^2 ( cos ( fc )4 c fc + 2Am s (£)4 +7r c fc ) , 

k— — 7C 

(31) 

where the mean-field staggered magnetization m s {t) has to 
be determined self-consistently. For developing an intuition 
it is worthwhile to imagine the dynamics of pseudo-spins 
in k-space by defining pseudo-spin operators of = ctcfc — 

cl +7l Ck+TT and a'l = c{ +7r c k + c\c k+n . Note that these 
momentum-space pseudo-spins are different from the orig- 
inal spins on the chain. In pseudo-spin representation the 
staggered magnetization is given by the average x-projection 

per pseudo-spin, m s = Y^h=-n/2 ( CT fc)' an( ^ tne mean -field 
Hamiltonian can be written as 

tt/2 

H MF {t) = -J J2 (cos(k)a z k + 2Am s (t)at), (32) 

fc=-7r/2 

The Neel state as an initial condition corresponds to all 
pseudo-spins pointing in x-direction at t = 0. Then they start 
to precess due to a Zeeman field that depends on the instan- 
taneous average orientation of the ^-projection of the spins. 
In these terms it is easy to understand the evolution of the 
staggered magnetization m s (t) for A = 0. We simply have 
a collection of independent pseudo-spins subject to constant 



Zeeman fields J cos k in the z-direction. Because the field 
magnitude varies from spin to spin over a bandwidth, they 
precess at frequency 2 J. Since the band of precession fre- 
quencies is continuous, the spins gradually dephase leading to 
the l/yi decay (15) of the oscillation envelope of m s (t). 

The situation is more complicated in the case ol A / 0. 
Now there is also a field in the a>direction, which is the same 
for all spins but changes in time according to the instanta- 
neous orientation of the spins. To lowest order in A, i.e. set- 
ting m s (t) = mo(t) = Jo(2Ji)/2 in the Hamiltonian (32), 
the additional Zeeman field in x-direction tilts the precession 
axis, giving rise to a smaller average x-projection of the spins 
and thus leading to a faster decay of m(t). The numerical 
results for the time evolution of the staggered magnetization 
according to (32) are shown in Fig. 10 for different values of 
A. Finite A leads to accelerated dephasing of the oscillations 
very much like in the unbiased calculations (Sec. III). How- 
ever, the asymptotic law as extracted from the numerical so- 
lution by Hastings and Levitov [97] for < | A| < 1 exhibits 

2 

algebraically decaying oscillations with aTs envelope, 

m s (t) ~ {d cos(wii + 4>i) + C 2 cos(w 2 i + 2 )} ,(33) 

Wi = 2 J, oj 2 = vl — A 2 . This algebraic decay, as well as the 
two frequencies, which lead to a revival phenomenon [97], is 
in contradiction with the MPS calculations for the full Hamil- 
tonian (1). For A > 1, the staggered magnetization saturates 
to a nonzero value for t — + oo, which is presumably also an ar- 
tifact of the mean-field approach not corroborated in the unbi- 
ased treatment. We conclude that the approach provides only 
a very rough picture of the order-parameter dynamics, which 
confirms the importance of collective effects, apparently not 
captured by the effective non-interacting mean-field Hamilto- 
nian (31). 
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V. XZ-MODEL - EFFECTIVE DESCRIPTION FOR A > 1 



In this section we study the time evolution of the staggered 
magnetization m s (t) following a quench from the Neel state 
in the analytically treatable XZ model. This serves as a com- 
plementary analytical approach to the numerical investigation 
of the quench dynamics in the XXZ model in the regime of 
large anisotropics A 3> 1 and allows of a discussion of the 
long-time asymptotic behavior of m s (t). The XZ model is 
defined by the Hamiltonian 



h xz = ./^ {2.s;.v;. , • a.s;.s;..} (34) 

3 

= H X xz + ^ X! {SfSf+i + S j S j+i} ■ 



At equilibrium the XZ model exhibits a quantum phase transi- 
tion at A = A c = 2 which separates two gapped phases with 
antiferromagnetically ordered ground states in z-direction for 
A > A c and in x-direction for A < A c . It differs from 
the XXZ model (1) by terms violating the conservation of 
Sfo t — S?, but has the advantage of being analytically di- 
agonalizable. In the following we will prove that the staggered 
magnetization in this model vanishes for all finite A > A c in 
the long-time limit after a quench from the Neel state and cal- 
culate the exact time evolution of m s (t) semi-analytically up 
to times Jt sw 100, thus going beyond the time window acces- 
sible by the MPS calculation for the XXZ chain. 

Using the Jordan-Wigner transformation for 5J and 5| and 
going over to momentum representation, the Hamiltonian of 




FIG. 11: Dynamics of the staggered magnetization m s (t) in the XZ 
chain initialized in a Neel state. Symbols correspond to numerical 
results, lines represent analytical results or fits by corresponding laws 
(see text). 



the XZ model (34) takes the form 
J 

k— — 7T 



tfxz = ! £ {(A + 2)cos(fc)4 



"A- 



+ -(A - 2) sin(fc) (alal k + a k a- k ) j , (35) 

with aJ k and a k denoting respectively creation and annihila- 
tion operators of spinless Jordan-Wigner fermions with quasi- 
momentum k. This Hamiltonian can be diagonalized by the 
Bogoliubov transformation, 



a-k 



cos 6u —ismtik 
-ismOk cos 8k 



a-k 



M A 



a-k 



2- A 

tan 26 h = tan k , 

2 + A 

which, maps (35) to a model of free fermions, 



(36) 



(37) 



with a dispersion Sk = J-\/l + A 2 /4 + A cos 2k. 

Since the initial Neel state is the ground state of (35) with 
A = A — > oo, it is convenient to express the time-dependent 
(Heisenberg) Jordan-Wigner fermion operators (t) in terms 
of the Bogoliubov quasiparticle operators a° which diagonal- 
ize the Hamiltonian (35) for the initial value A = Aq, 



a- k (t) 

4(*) 



Ma 









Ae k t 



M7 l M Ao 



(38) 



This reduces the computation of correlation functions at arbi- 
trary time to the evaluation of ground-state expectation values. 

To calculate the time evolution of the staggered magnetiza- 
tion m s (t) following a quench in the XZ model, we define the 
two-spin correlation function, 



(39) 




FIG. 12: Relaxation times and oscillation period T — — as a func- 
tion of anisotropy in the XZ model. Algebraic laws are emphasized 
by solid lines. 
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from which the square of the staggered magnetization is ob- 
tained by taking the infinite-range limit, 

m 2 s (t) = lim C{£,t) . (40) 

l— >oo 

In the fermionized picture of the XZ model the two-spin cor- 
relator takes the form 

(S 2 Sf) = ^(-l) e (A B 1 A 1 ...Bt_ 1 A^ 1 B e ) , (41) 

with Aj = ajj + aj and Bj = aj — aj being the Majorana 
operators at lattice site j [102]. Using Wick's theorem, this 
correlation function can be expressed as a Pfaffian of pairwise 
contractions [42]. For the quench problem studied here, the 
explicit form of these contractions follows from (38) and is 
given by 



(AjAi) = (BjBi 



(42) 



-ik(J—i) 



sin 2skt sin 2(f>k , for i ^ j, 



(AjBi) = r 

J — 7T 



dk 
2^ 

^ e -*a-i) e *ae fc(cos2 ^ 



2tt 

— isin20fe cos2£fct) , (43) 

with <f>k = Ok — 6® (see also [62], where identical expressions 
have been derived for the transverse-field Ising model). We 
have taken the thermodynamic limit and converted the sums 
into integrals in the expresions above. In the limit t — > oo 
for A > A c the evaluation of (40) reduces to the computation 
of a Toeplitz determinant, since the contractions (43) of the 
Aj's and Bj's among themselves vanish. Szego's theorem 
can then be used to calculate the asymptotics of the Toeplitz 
determinant, yielding the result 



lim C(£,t) 

t— >oo 



«>i 1/1 + V 1 - 4 / A 2 



(44) 



Thus, after a quench from the Neel state in the XZ model, 
the staggered magnetization vanishes for all finite A > A c at 
large times. 

At finite times, when the contractions (43) do not vanish, 
the Pfaffian representing the two-spin correlator (41) can be 
evaluated numerically at arbitrary times for a given distance. 
Due to the so called light-cone effect [53, 91], two spins 
at a distance I are not causally connected at times smaller 
than ut < i/2, since the correlation length of the initial 
Neel state is zero. Here u denotes the maximum (classical) 
speed of quasiparticles, which in the XZ model is given by 
u = maxfc(9fc£fc) = 2J. Exploiting this light-cone effect, 
the staggered magnetization can be calculated in terms of a 
finite-range correlation function, 



!(*) 



2Jt<4 



C(£,t) 



(45) 



This method significiantly reduces the computational effort at 
short times. We remark however that the light cone is not com- 
pletely sharp in quantum-mechanical systems [53]. Never- 
theless, for practical finite-precision calculations the infinite- 
range limit of the two-spin correlator is reached for distances 
just a few lattice sites beyond the light cone. 



The results for the time evolution of the staggered mag- 
netization following a quench from the Neel state in the XZ 
model are displayed in Fig. 11. As is the case for the XXZ 
chain, an explicit analytical expression for m s {t) in the XZ 
model can be derived for a quench to A = 0, which is given 
by m s (t) = 0.5 cos 2 (Ji). For A < A c , the numerical data 
for m s (t) at large times fits very well exponentially decaying 
oscillations of the form 



m s {t) oc e */ T2 (cos 2 (u)t) — const.) 



(46) 



In this regime, the behavior of m s (t) in the XZ model is qual- 
itatively different from that in the XXZ model, as can be seen 
from the period of the magnetization oscillations. In the XZ 
model the period diverges at the critical point (see Fig. 12) 
, whereas it becomes smaller upon approaching the isotropic 
point in the XXZ model (see Fig. 4). Furthermore, the criti- 
cal point exactly marks the crossover between oscillatory and 
non-oscillatory behavior of m s (t) in the XZ model. 

For A > A c , the staggered magnetization decays exponen- 
tially in the XZ model and shows no oscillations at large times. 
Interestingly, the numerical results for m s (t) in the XXZ and 
XZ models are almost indistinguishable at large anisotropies 
A ^> 1, as can be seen from Fig. 11. We have extracted the 
relaxation times from exponential fits to the numerical data, 
obtaining a clearly pronounced minimum right at the isotropic 
point (see Fig. 12). The relaxation time scales as T2 oc A -1 
for A < A c and as n oc A 2 for A > A c . 



VI. GAPLESS THEORY - LUTTINGER MODEL 

In the analysis of the XX limit (Section II) it became ev- 
ident that, if the initial gap is sufficiently small, the non- 
oscillatory relaxation of the order-parameter dynamics is de- 
termined by low-energy modes, which motivates the applica- 
tion of the Luttinger model to a quench to the gapless phase 
A < 1 of the XXZ model. 

In Section A the Luttinger model, -Hll = 
£ J dx {if (TiTIfz)) 2 + ± (d x <j>{x)f\, has been intro- 
duced as a low-energy effective theory for the XXZ chain in 
the easy-plane regime. The bosonized form of the staggered 
magnetization is given by m s ~ (cos(2^>)) ;r= o, where made 
use of the translational invariance. The remaining problem 
amounts to computing the time evolution of (cos (20)), 
starting from a state where the field <fr is initially pinned at 
or tt/2. We remark that this problem is essentially the 
dual of the dephasing problem studied in [103], and thus 
we expect an exponential decay of m s with a characteristic 
time scale r ~ J/(KA S ). A convenient technique for 
solving this problem is the truncated Wigner method [59], 
which is exact for quadratic Hamiltonians such as (B9). 
Using this approach, the time-dependent expectation value of 
the staggered magnetization can be written as a functional 
integral over the Wigner transform Qw{4>o, 4>o) of the initial 
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density matrix: 



VII. GAPPED THEORY - THE SINE-GORDON MODEL 



(cos(20)) = / V<f>(t) J V(ct>o,4>o)Qw((t>o,<t>o) 

x cos(20) 5(<f> - u 2 dl<j}) (47) 
T>((f> ,4>o) Qw(4>o, <t>o) cos(20 d (x, tj) 



Here, the functional (^-distribution ensures that one integrates 
only over solutions of the classical equations of motion and 
<f> c i(x,t) denotes the classical solution of the ID wave equa- 
tion corresponding to the initial conditions <fio(x) and cbo(x). 
We have also used the fact that the operator cos(</>) is diagonal 
in the (^-representation. The solution (f) c i(t) can be explicitly 
constructed using d'Alembert's formula. After switching to 
dual-field representation using Kud x 6 = </>, we get 



(cos(20)) ~ / V((f> ,9 ) g w (0 O A) cos |0 O (x - ut) 

+ (j) (x + ut) + K9 (x + ut) - K6 (x - ut) 



Since in the initial state <j> is pinned at (f> = 0, we factor out 
the c/> dependent part of the integral, obtaining 



m s (t) ~ (cosK(8(ut) - 9(-ut))) , 



(48) 



where the brackets with the index denote the expectation 
value taken with respect to the initial state. The r.h.s. of Eq. 
(48) can be estimated within a semiclassical analysis, where 
the ground state of the Luttinger Hamiltonian (B9) with an 
additional mass-term ~ A s (j> 2 is used as the initial state. This 
finally leads to 



m s (t) 



cxp-^-((6(ut)-e(-ut)) 2 ) c 



cxp —K I dq 



■ sin 2 (qui) 



A s t» 



exp(-wKA s t/2) , 



(49) 



where A s again denotes the gap of the initial state[177]. In 
contrast with the empirical rule (30) for the XXZ model, the 
Luttinger model, being a continuum theory, does not repro- 
duce oscillations. The non-oscillatory relaxation in (49) is 
characterized by a relaxation time inversely proportional to 
the gap and to the Luttinger parameter, r 



tKA, 



a behavior 

identical to the conformal field theory result (8) and similarly 
observed in the numerical calculation for the quench in the 
XXZ model. However, the algebraic prefactor present in the 
case of the XXZ model (30) and the spin-density-wave ini- 
tial state under the XX Hamiltonian (26) is not present in this 
treatment of the Luttinger model. Since the Luttinger model 
includes the XX limit at K = 1, we conclude that the miss- 
ing algebraic prefactor is a shortcoming of the initial state, 
which has been approximated as the ground state of the Klein- 
Gordon model (B16). More accurate results could provide a 
treatment using the sine-Gordon Hamiltonian, which as we 
shall see in the next section strongly complicates the problem. 



In this section we analyze the quench in the sine-Gordon 
model (see appendix B), 



Hsg = ^ [ dx[uK(nIl(x)) 2 + ^(V<f>(x)) 2 ] 
. dcccos(40(.T)), 



2JA 

(27ra) 2 



(50) 



as a possible continuum approach to the quantum quench the 
XXZ chain for A > 1. In what follows we use the boundary- 
state formalism as a convenient tool for describing the non- 
equilibrium problem [55]. In this formalism the initial state, 
which is not the eigenstate of the Hamiltonian, can be thought 
of as a special superposition of pairs of eigenmodes of the 
quantum Hamiltonian with opposite momenta, which sums up 
into a squeezed state of eigenmodes [ 1 04] . Unlike for the gap- 
less Luttinger-liquid theory, we cannot present a full solution 
of the dynamics. Possible directions to be followed in future 
are pointed out. 

Since the sine-Gordon model has relativistic (Lorentz) in- 
variance, we can exchange the time and space directions x <-> 
t and consider the following boundary-in-time Hamiltonian 
(using more conventional notation x again for the imaginary 
time direction) 



H = 



1 

2^ 



dx[uK(TrU(x)f 



K 



(V</>(*)) 2 



2 J 



(27ra) 2 
A cos(2^(x))(5(x)} , 



dx{Acos(4(l>(x)))6(x) 



(51) 



where 8(x) is a theta-function and S(x) takes care of the ini- 
tial condition. In order to implement the Neel state as an ini- 
tial condition we send Ao — > oo which corresponds to the 
Dirichlet boundary (initial) condition. This boundary (initial) 
condition formulation can be reformulated in a boundary-itafe 
formalism of the boundary sine-Gordon model (bSG). The ini- 
tial condition is expressed as a squeezed state of bulk degrees 
of freedom. We note that for noninteracting particles or Lut- 
tinger liquid this correspondence can be seen directly. Since 
for K < 1/2 there are only solitons and antisolitons in the 
spectrum (repulsive regime of the sine-Gordon model), we ob- 
tain the boundary state in the following form 



\B(t = Q)) D = Afexp 



K${e)A\{e)A\{- 



|0).(52) 



Here A\ b (6) is an operator of creation of the soliton (a) or 
antisoliton (b) and Kf^(6) is a reflection matrix of soliton- 
antisoliton pair corresponding to the Dirichlet boundary con- 
dition. The rapidity 6 is related to the momentum P = 
M s sinh 9 and energy E = M s cosh 9, where the soliton mass 
M s is given by [105] 



r(£) , 



2r(|) 



v/HXl + f) 



(53) 
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where we define £ = /3 2 /(l — 1 ). 

The evolution is trivial in the soliton basis, because the bulk 
Hamiltonian is diagonal in soliton-antisoliton operators, 

\B(t>0)) D = Afexpfj K$(6,t)Al(6)Al(-e)) |0) 

K^(0,t) = K$(6) cxp(2itM s cosh(0)) . (54) 

In the boundary state formulation the evolution of the mag- 
netization is equivalent to the computation of the following 
quantity 

m s (t) = (B(t)\ cos(20(O))|B(i)> . (55) 

In general the squeezed state represented by the boundary 
state \B(t)) should be expanded as a series in powers of the 
reflection matrices. This produces multiple dynamical pro- 
cesses which include solitons and antisolitons. Multi-particle 
expectation values of the operators, like cos(2(/>), are called 
form-factors. To compute the correlation functions in the mas- 
sive theories at equilibrium, only a small number of lowest 
form-factor contributions is necessary. However, our evalu- 
ation of the lowest order contributions in our case provided 
results contradictory to the numerical simulations. The reason 
for this will be found in the spectral analysis of the next sec- 
tion, which hints that not only soliton-antisoliton form-factors 
are important (which is the case for the spectral function of 
the sine-Gordon model for small energies), but also multiple 
processes, which include energies well above the spectral gap 
(soliton mass), are necessary to be considered. Technically, 
the problem of inclusion of multi-soliton form-factors is rather 
difficult. The ^-integrals corresponding to evaluation of dif- 
ferent multi-particle contributions become even more compli- 
cated because of the reflection matrices K c j^(ff). 

A possible alternative approach to this form-factor evalua- 
tion could be a resummation of the leading divergencies of the 
scattering processes in the presence of a boundary state. Since 
the ultra-violet energies give an important contribution in our 
problem, one can try to proceed by considering the logarithm 
of the one- or (two-) point function and to sum the leading 
contributions as proposed previously [106-109]. However, 
the complexity of the boundary reflection matrix does not al- 
low to realize this program. We hope to return to this problem 
in future. 

In view of the high complexity of the boundary-state for- 
malism it may be worthwhile to establish phenomenological 
analogies between non-equilibrium dynamics and equilibrium 
dynamics of the sine-Gordon model. This would be useful 
since for calculating dynamical structure factors a powerful 
machinery has been developed over the last decades. Argu- 
ing that the initial state can be described by a thermal ensem- 
ble (with some effective temperature considered as a fitting 
parameter) instead of the boundary state, we can relate the 
dynamics of the magnetic order parameter to the two-point 
function, 

m s {t) ~ (cos(20(t)) cos(20(O))) , (56) 

where the average is taken over some thermal ensemble char- 
acterized by temperature T e ff. The operation of cos(20(O)) 



onto the thermal state is a possibility to introduce some mag- 
netic order, or, stated otherwhise, to establish an analogy to 
Eq. (55); cos(20(t)) acting on the thermal state is a way 
to mimic a boundary-in-time state. We note that such ap- 
proach has been successfully applied for studying dynamics of 
a non-local observable in quench in the quantum Ising chain 
[110]. The dynamics of the two-point function (56) is sepa- 
rated into two regimes: large-temperatures T e tf M s and 
low-temperatures T e ff <C M. It is known that for large 
energies (UV) massive models, like the sine-Gordon model, 
have the conformal filed theory asymptotics. Therefore, in the 
large-temperature regime the behavior of the correlation func- 
tions should be the same as in the high-temperature limit of the 
corresponding conformal field theory. Hence, for T e f f ^S> M s 
the large-time asymptotics of the correlation function is given 
by an exponential decay 

m s (t) ~ cxp[-7rT e// yi]. (57) 

This conformal field theory behavior is universal also for the 
gapless phase, where, at least in some regimes of weakly mag- 
netized initial states, setting T e ff = A s this behavior is a 
good first approximation of the dynamics of the order param- 
eter in the quench problem (see Table I). However, in the 
gapped phase we cannot find a reasonable way to define T e ff . 
For example, the temperatures corresponding to the Boltz- 
mann ensembles used in the following section do not repro- 
duce at all the numerical findings. 

In the other regime, T e fj <C M s , the structure of the mas- 
sive theory is important. In this case the leading order behav- 
ior comes from the zero-momentum exchange processes and 
depends on the structure of scattering matrix £(0) in this limit. 
Resummation of the kinematical singularities leads again to 
the exponential decay for the two-point correlation function 
[111], in agreement with a quasi-classical formula from Ref. 
[112]. Implementing results of [111] to our situation we ob- 
tain 

m s ~ exp[-T e f f e- M '/ T 'tn], (58) 

where the proportionality coefficient depends on the power 
of M s . Such behavior however is in disagreement with our 
numerical findings, where in the limit of large M s (large A) 
we find a decay rate proportional to A~ 2 . 

We conclude that although the sine-Gordon is a valuable 
candidate for describing the dynamics following a quantum 
quench in the XXZ model, the evaluation of the correspond- 
ing form-factors is difficult and demands further efforts. A 
relation of the coherent dynamics of the order parameter to dy- 
namical structure factors, circumventing this problem, is not 
straightforward to be established. 

VIII. SPECTRAL ANALYSIS 

For a deeper understanding of the relaxation dynamics, it is 
useful to consider the problem in energy space. The idea is to 
associate properties of the spectrum of the Hamiltonian to the 
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FIG. 13: Analysis of the spectrum of the XXZ chain for a system of 14 sites with periodic boundary conditions. Dotted lines mark the position 
of the energy of the system, E = (il>o\H |V>o). See text for the description of the regions marked (i)-(iv). 
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FIG. 14: The frequency distribution of the staggered magnetization, /m„(e), for a system of 14 sites with periodic boundary conditions. 
Histograms resolve individual peaks, solid lines correspond to (non-unique) smoothened distributions emphasizing separable contributions. 
For A = the exact solution in the thermodynamic limit (14) is drawn instead. 



dynamical phenomena observed in the simulation of the time 
evolution and to clarify the possibility of separating energy 
scales - a question which is especially important for improv- 
ing analytical descriptions of the non-equilibrium dynamics. 

Using the Lehmann representation, the time evolution of 
an operator O takes the form of a Fourier transform over the 
eigenlevels of the Hamiltonian, 



which, via the distribution function 



<o(*)> = E' 



-it(E m -E„ 



<(i/j \m)(m\O\n)(n\i/j ) 



For a more convenient continuum description, we introduce 
the quenched probability distribution, 



PVoO) = J2 5 ( e + E o- e n)\{n\ipo}\' 



(59) 



which determines the properties of the stationary state at t — > 
oo [78, 113-115] (the frequencies e are shifted by Eq - the 
ground-state energy of H). It can be compared to the thermal 
(Boltzmann) distribution of the grand canonical ensemble, 



-y 



S(e + E - e n )e 



,/T 



(60) 



where the temperature is set by the energy of the initial state, 
/ deepsie) — (i^q\H\i[;q). In general, it is known that the 
thermal distribution can deviate strongly from the quenched 
distribution [63, 78, 113-1 15], which leads to the phenomena 
of absence of thermalization. Here the thermal distributions 
are used only as a reference and questions in connection with 
thermalization phenomena will not be investigated. While the 
quenched probability distribution, p^ , captures the effect of 
the initial state, the distribution of the expectation value, 

0(e', e)=J2 6 ( e ' + E o- tm)5(e + E - e n ){m\0\n) , 

reflects the specific spectral properties of the given observable. 
p^ and 0(e', e) provide the contributions to the weighted ex- 
pectation value, 

W (e', e)=J2 S(e'+ E Q - e m )6(e + E Q - e n ) 

n,m 

x(-0oM (m\0\n)(n\ipo) , (61) 



fo(e) = J de'Wo(e + e',e'), (62) 
represents the dynamics of the observable in frequency space, 
0(t)= / dee- lit fo(e). (63) 



The spectral properties of the XXZ chain prepared in the 
Neel state, |^>o) = |f J,f • • • ||), with the staggered magneti- 
zation as the observable, O = m s , are calculated by means 
of exact diagonalizations for small system sizes. Fig. 13 dis- 
plays the results for a chain of length = 14. The small sys- 
tem size results in strongly peaked distributions, but we have 
made sure that the qualitative features we extract in the follow- 
ing analysis are stable against variations of the system size, 
both towards larger, N = 18, and smaller, N = 10, values. 
Quantitative information cannot be extracted from this simple 
analysis, but this might be possible when going to larger sys- 
tem sizes by means of more involved techniques, such as the 
Lanczos method [71]. 

In the non-interacting limit (A = 0), where the Hamilto- 
nian has a free-fermion representation as discussed in section 
II, all distributions are centro-symmetric about e = Eq- The 
quenched distribution (e) exhibits peaks at e — Eq = ± J, 
which are not present in the thermal distribution. From the 
discussion in section II and the finite size-study [Fig. 13(a-c)] 
we can separate two contributions to W ms (e', e): 

(i) Free-fermion band edges. The distribution has maxima 
and sharp cutoffs at the edges e — e' = ±2 J. 

(ii) Low frequency. The contributions from the areas at e — 
e' <C J are weaker than the band-edge contributions 
(i). 

For small but finite anisotropies, as exemplified by the results 
for A = 0.4 [Fig. 13(d-f)], the main features of the A = 
distributions are still present. However, the distributions lose 
their reflection symmetry and the weight is shifted towards 
lower energies. In addition to this asymmetry, the cutoff of 
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the distribution of m s is no longer absolutely sharp at the band 
edge. This can be understood as an effect of the breakdown of 
the free-fermion quasi particle description of the Hamiltonian 
at A ^ 0. 

For large anisotropies in the gapped regime [e.g. A = 3, 
Fig. 13(m-o)], the quenched probability distributions have a 
single peak at e = and a continuum above the gap, which 
is rather fiat in comparison with the thermal distribution. The 
expectation values of the observable m s (e', e) are spread over 
a large energy region, however, the relevant weighted expec- 
tation values W ms (e',e) receive contributions of only three 
different types, all of which are located at low energies: 

(iii) Continuum above the gap. The contributions form 
(m\m g \n), n ^ 0, to ^ form a rather continuous 
distribution with frequencies far below the gap e — e' <C 
JA. 

(iv) Off-diagonal contributions involving the ground state. 
Elements (0\m s \n), n / 0, result in frequencies equal 
to or larger than the gap. 

(v) Ground state. An isolated peak is located at the ground 
state energy of H (e = e' = 0). 

For lower A this picture remains basically valid. However, 
in Fig. 13(j-l), where the results for A = 2 are shown, the 
width of the contribution (iv) becomes of the order of the gap. 
For even lower anisotropies, e.g. A = 1, separation of the 
contributions is no longer possible [Fig. 13(g-i)]. 

The effects of these contributions on the frequency distri- 
bution fm B ( e ) i s shown in Fig. 14. The characteristics of 
time evolution can now be inferred from the properties of 
the Fourier transform of the distribution f„ ls (e). The alge- 
braically decaying oscillations at A = (13) are a conse- 
quence of the step-like shape of /m 3 (e), which is a conse- 
quence of the properties of the contribution (i). For A > 
the smearing of the edge is the cause of the exponential de- 
cay of the oscillations. The finite width of the edge in Fig. 
14(b) is consistent with the relaxation time of the oscillations 
(Fig. 4). At A > 1 there are also contributions at large fre- 
quencies originating from contributions of type (iv), resulting 
in rather broad peaks in f ms (e), which lead to the quickly 
decaying oscillations seen in Fig. 5. The low-frequency con- 
tributions (iii) are reflected in a peak of f ms (e) around e = 0, 
whose width corresponds to the relaxation time n of the non- 
oscillatory decay (Fig. 4). 

The isolated peak (v) at zero energy is irrelevant for the 
dynamics. It would correspond to a finite asymptotic value 
at t — > oo, which apparently vanishes in the thermodynamic 
limit. 

In summary, on the basis of the analysis of the frequency 
distributions, we can now draw a qualitative crossover pic- 
ture from the oscillatory to the non-oscillatory behavior as A 
varies. Approaching the isotropic point from small values of 
A, the band edges (i) are smeared out, leading to a decreas- 
ing relaxation time T2. Starting from large values of A, the 
low-frequency peak merges into the homogeneous distribu- 
tion upon approaching the isotropic point. Hence the char- 
acteristics of (i) and (iii) contributions, dominating at small 



(respectively large) values of A, are both lost at intermediate 
values of A, where the interplay of all energy scales appar- 
ently leads to a non-generic dynamical relaxation of the order 
parameter. We also note that, even in the regime A 3> 1, 
where the initial state is rather close to the ground state of the 
Hamiltonian, the relevant part of the spectrum located above 
the gap is a multi-particle continuum, difficult to be treated 
analytically. 



IX. CONCLUSIONS AND OUTLOOK 

We have analyzed the dynamics of the staggered magneti- 
zation in quantum spin chains following a quantum quench, 
considering various antiferomagnetically ordered initial states 
by using a number of complementary numerical and analyti- 
cal approaches. In the numerical MPS study we have essen- 
tially found three types of relaxation dynamics for the order 
parameter: (i) For highly ordered initial states (Ao ^> 1) and 
sufficiently small anisotropy parameters of the Hamiltonian at 
t > there are Rabi oscillations, which dephase exponentially 
in time away from the XX limit; (ii) for strong anisotropies 
(A ^> 1) there is an exponential decay and the relaxation time 
scales as r oc A 2 ; (iii) for initial states close to the phase 
transition we found evidence for algebraic corrections to the 
exponential decay. There is a crossover phenomenon between 
oscillatory (small A) and non-oscillatory dynamics (large A), 
but no clear point of transition can be identified. Either both 
types of dynamics superimpose on each other, as it is the case 
for Ao > 1, or both vanish in an extended transition regime, 
in which a non-generic behaviour is found (case of Ao ^> 1). 

We have shown that a precise description of the Rabi os- 
cillations (i) is possible only when the full spectrum is taken 
into account. Therefore mean-field as well as low-energy 
approaches lead to incomplete results - an analytical treat- 
ment of this type of dynamics is feasible only by novel ap- 
proaches. It has become clear that quasi-particles relevant for 
Rabi-oscillations are located at the band edges. In contrast to 
equilibrium properties, where many -body effects can be incor- 
porated into the relevant linearizable Fermi-level excitations 
within the bosonization formalism, it is not obvious how to 
treat interactions in combination with the quadratic dispersion 
relation at the edges of the band. A treatment along the lines 
of previously developed concepts dealing with non-linear ef- 
fects in dynamical phenomena [50, 116] might be a possible 
solution to this problem. 

While the exponential decay at large anisotropies (ii) ap- 
pears to be a rather generic behaviour and is also reproduced 
in the exactly solvable XZ model, the algebraic prefactor 
in front of the exponential law (iii) is a more intricate phe- 
nomenon. In the XX limit, where we approximated the initial 
state by a spin-density wave, such order-parameter dynam- 
ics are reproduced. However, standard field-theoretical ap- 
proaches, such as conformal field theory or the description by 
Luttinger model adopted here, capture roughly the scaling of 
the corresponding relaxation time, but do miss the prefactor. 
Our results suggest that this is an consequence of the phe- 
nomenological description of the initial state in terms of a sim- 
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pie massive theory (Klein-Gordon). A more elaborate treat- 
ment of the initial state using the sine-Gordon model could 
resolve this deficiency of the field-theoretical descriptions. 

The sine-Gordon as well as the original XXZ model is in- 
tegrable. Analyzing the quench by using the integrability of 
these models amounts to evaluation of form factors of parti- 
cles with non-trivial statistics. While for two or three particles 
this problem can be solved [117, 118], we have found that 
non-equilibrium dynamics require the evaluation of higher or- 
der form factors - a yet unsolved and highly complex prob- 
lem. A promising approach is to use the structure of the 
Bethe-ansatz solution in combination with a numerical algo- 
rithm [67]. 

Similarly to studies of other models [84, 119] we find as 
a generic feature that relaxation times become small in the 
vicinity of the quantum critical point. However, the behaviour 
is far too rich to be attributed to a generic dynamical phase 
transition [119]. In our analysis we have shown that effective 
descriptions by low-energy theories belonging to the univer- 
sality class of the model can not capture all relevant processes. 
A sort of a dynamical phase transition occurs, however, in the 
mean-field description (a finite magnetization is found in the 
long-time limit for A > 1), which treats interaction terms as 
on an infinite-dimensional lattice. This hints that the existence 
of a sort of dynamical critical behavior may, very much like 
for equilibrium phase transitions, depend on dimensionality. 
A first step towards understanding of the role of dimension- 
ality is the study of non-equilibrium dynamics in infinite di- 
mension, for example by using dynamical mean-field theory 
[86, 119]. How to treat coherent dynamics in a two- or three- 
dimensional system is still an open question. 

Finally, we have to mention that the Heisenberg chain is a 
simplified model, well suited for numerical and analytical in- 
vestigations, but not necessarily appropriate for full descrip- 
tion of experimental systems. Although in experiments with 
two-level atoms in optical lattices behaviour similar to our 
theoretical prediction is observed [19, 21], the model has to 
be adjusted to provide an accurate description of the experi- 
mental results. For example, the effect of density fluctuations 
beyond the purely magnetic model needs to be investigated. 
Although the matrix product algorithm is an efficient method 
to study the relaxation dynamics, the numerically reachable 
times are fundamentally restricted by growing entanglement. 
The runaway time may become very small in models more 
involved than spin-^ chains, in particular when particle fluc- 
tuations need to be taken into account [80]. Recently, schemes 
have been proposed which allow to go beyond what is possi- 
ble within the conventional matrix product algorithms [ 1 20- 
123]. The main idea in these approaches is to calculate di- 
rectly the dynamics of an observable rather than explicitly 
follow the evolution of the state. Another important aspect 
neglected here, but relevant in experiments, is temperature. 
How to efficiently include effects of finite temperature within 
a time-dependent matrix product algorithm is a yet unsolved 
problem [121, 124]. 
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APPENDIX A: QUANTUM MAGNETISM IN OPTICAL 
LATTICES 

The underlying model for realizing quantum magnetism in 
optical lattices is a single-band Hubbard model, 

H = {''..•" + H - c -} + U n n H n il 

ija i 

+ J~] -7r-(?W - - y^Mio-nio- . (Al) 

ia i 

a\ a are the creation operators of a particle in a Wannier state 
of type a at site i, satisfying bosonic commutation relations. 
nia = a \ a ai<j are the corresponding occupation numbers. The 
internal degree of freedom a represents typically a hyperfine 
state and can be identified with (pseudo) spin-i, a = f, j. 
The hopping integral i; Jcr and the interaction parameters U aa i 
depend on the geometry and the depth of the lattice and can be 
expressed in terms of overlaps of the Wannier orbitals. If, for 
concreteness, we consider a periodic, spin-dependent lattice 
potential with an isotropic spacing a, 

V a {v)= J2 V^wpKr^, (A2) 

with K = — , the recoil energy needs to be much smaller than 

the lattice depth, E r = 2 ^ -C V^ a , so that the atoms are in 
the lowest harmonic level and single-band description (Al) is 
valid. In reality, the gaussian shape of the laser beams intro- 
duces inhomogeneity in the lattice depth in addition to the har- 
monic trapping potential The superposition of polarized 
laser beams generates spin-dependent potentials [31, 126]. In 
general, the orbitals extend only over short distances and the 
hopping can be restricted to nearest neighbors, = t^ a , if 
rj — v,j = e^a, where [127] 

V * (4/^)S ) L / 4 (V) 3/4 exp[-2(V/^) 1 / 2 ](A3) 

U aa , « (S/^/^fcawX^TwTV^VW) 1 ^) 

Here, V^- a = iV^V„J (V^ 2 _+ V*( 2 ) 2 , is the spin- 
average potential in each direction, V Mcrcr = V^, and a sacr > is 
the s-wave scattering length between atoms of spin a and a' . 
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In the case of strong on-site repulsion, t^ a <C U aa ', and 
integer filling, the system is in the Mott phase, where occu- 
pation number fluctuations are essentially suppressed. In this 
case, the effective basis contains locally only singly occupied 
spin-up |t) or -down |J.) states (for the sake of simplicity we 
choose + = 1, although higher occupation num- 
bers are also possible). In this subspace, neglecting terms of 
order tfj/U^ a ,, the Hubbard model (Al) can be mapped onto 
a spin-i Heisenberg model (XXZ model) [128-131], 
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FIG. 15: The ground-state phase diagram of the XXZ model in one 
dimension. 



The superexchange interaction constants J l J and Jj 5 are given 
by 
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(A6) 
(A7) 



An analogous treatment can be carried out for the fermionic 
Hubbard model. In the resulting magnetic Hamiltonian (A5), 
J± has the opposite sign compared to Eq. (A7) and in the 
expression for J z the last two terms are absent since double 
occupancy is forbidden by the Fermi statistics. 

For appropriately chosen lattice and interaction parameters, 
the anisotropy of the spin exchange, A 1 - 7 = Jj 7 / is tun- 
able to a large extent. For example, in the bosonic case with 
symmetric on-site repulsions, a ferromagnet with possible 
easy-axis anisotropy, A' y = ^(t^- + t 1 ^) > 1, is realized. 
In addition, as demonstrated recently [19], double- well poten- 
tials can be used to change the sign of the exchange interac- 
tions from ferro- (,Pj < 0) to antiferromagnetic { J l J > 0) 
[80]. 

Although the Heisenberg Hamiltonian (A5) is a good first 
approximation to strongly interacting two-component Bose- 
gases in optical lattices, we note that the measurements of 
Trotzky et al. [19] clearly show the limitations of the purely 
magnetic picture. The strong repulsion leads to superex- 
change interaction which reaches the order of currently re- 
alistic temperatures, J/fcs ~ 10 -9 -ftT, and in non-equilibrium 
experiments the dynamics slow down, so that effects of inho- 
mogeneous laser beams become strong. For larger tunnelings 
density fluctuations are important and introduce an additional 
higher frequency; excitations to higher Bloch-bands may also 
become possible. We conclude that, although the experimen- 
tal progress looks promising, further improvements in exper- 
imental setups are still needed in order to produce a clean re- 
alization of a quantum magnet. 



correlated systems in general, can be realized experimen- 
tally as proposed in the preceding Section. Here we give an 
overview of the equilibrium phases of the Heisenberg chain, 
focusing on antiferromagnetic exchange interactions. At the 
same time, the important concepts and notations to be used 
in the ensuing discussion of the non-equilibrium problem are 
introduced. 



The one-dimensional spin-i Heisenberg chain, 



(Bl) 



is integrable - the eigenstates and an infinite number of con- 
served operators can be obtained using the Bethe ansatz [132— 
136]. A number of equilibrium properties can be exactly com- 
puted for the Bethe wave function - examples are the energy 
and momentum of low-lying states, or local observables such 
as the staggered magnetization [137, 138]. For some specific 
cases, non-local properties can also be calculated analytically 
[139, 140] or by means of a combination of the Bethe ansatz 
with numerical algorithms [67, 141, 142]. A simplified in- 
sight into the physics of the Heisenberg chain can be gained 
from a continuum description via the bosonization technique 
[101, 143]. Here, results from both approaches, Bethe ansatz 
and bosonization, will be presented. 

The ground-state phase diagram of the XXZ model is rep- 
resented in Fig. 15. Without loss of generality the coupling 
J can be considered to be positive and the phases are simply 
characterized by A. The long-range ordered antiferromag- 
netic phase for Ising-like anisotropies A > 1 exhibits a spec- 
tral gap. In the easy-plane regime |A| < 1, a critical gapless 
phase is found. The phase for A < —1 is ferromagnetically 
ordered. 

A useful equivalent representation of (1) is a model of in- 
teracting spinless fermions, 



H 



J 



XXZ = 



E{4 



c ]+l c j 



:Act Ci ct +1 c j+1 }(B2) 



APPENDIX B: EQUILIBRIUM PROPERTIES OF THE XXZ 
MODEL IN ONE DIMENSION 



In a spatially anisotropic optical lattice the Heisenberg 
chain, a paradigm in the theory of magnetism and strongly 



obtained from (1) by Jordan- Wigner transformation from spin 
operators to spinless fermion operators [144], 
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(B3) 
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In the case of a one-dimensional Hamiltonian with nearest- 
neighbor interactions, particle statistics is irrelevant, and alter- 
natively the fermions can also be replaced by hardcore bosons 
[101]. 

The fermionic picture is especially useful in the non- 
interacting case (A = 0, also known as the XX limit), where 
(B2) is diagonal in Fourier space, 



Hxx = ^ £fec;.c fc , 
k 

6k = — Jcosfc. 



(B4) 



In the case of zero magnetization, which is of interest here, 
the ground state is described by the half-filled Fermi sea, 



xx 



1 10> , 



n ^ 

-7r/2</c<7r/2 

where |0) is the fermionic vacuum, Cfc|0) 
dinal) spin-spin correlation function, 

1 E 



(B5) 
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0. The (longitu- 
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which characterizes magnetic ordering, can be calculated ex- 
actly in the XX limit [102, 145-150]. The result is a super- 
position of quasi-long-range ferromagnetic and antiferromag- 
netic correlations, decaying by a power law, 

i-i-iy 



G zz {£) oc 
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(B7) 



For finite A, the extraction of correlation functions from the 
Bethe ansatz solution is highly non-trivial and only possible 
for some special cases (e.g. [140]). 

In order to obtain a continuum description of the Heisen- 
berg chain, the spectrum of the non-interacting model is lin- 
earized at the Fermi points and the modes are separated into 
left- and right-movers, 
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The cutoff A is of the order of the bandwidth. Starting 
from (B8), interactions can be included using the bosoniza- 
tion formalism [101]. At the renormalization-group fixed 
point, which captures the long-distance properties, the Lut- 
tinger model, 



H 



LL 
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dx{ K(irIL(x)Y 



— {d x (j){x)f 



(B9) 



provides the effective description for | A| < 1. IT(.t) and 4>(x) 
are conjugate bosonic fields, [n(x), 4>(x')] = iS(x — x'). We 
note that the excitations described by the Luttinger model cor- 
respond to linearly dispersed spin waves with velocity u. The 
values of both, u and the Luttinger liquid parameter K, can be 
derived from the Bethe ansatz [151], 



K = 



(B10) 



1 

Jsin(7r(l -/3 2 )) 
2(1-/32) ■■ 



where (3 is determined from the relation A = — cos 7r/3 2 . The 
functions AT (A) and it (A) in the antiferromagne tic regime are 
plotted in Fig. 16, K(0) = 1, K(l) = |, u(0) = J and 
m(1) = ^, additionally AT(-A) = K~\A). The bosonic 
fields can be mapped back to the spin operators, 

1 (-D x 

S z (x) = V4>(x) + ± — ^cos(20(x)), (Bll) 

7r ira 



where a ~ 4-. Here, the lattice spacing is set to one, so that 
the original sites are located &tx = i,i = l,...,N (N being 
the number of lattice sites). 

For the quadratic Luttinger Hamiltonian (B9), the correla- 
tion functions can be evaluated [101], 
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FIG. 16: The Luttinger parameter K, the velocity u, the gap A s , 
and the staggered magnetization m s in the XXZ model, calculated 
by Bethe ansatz. 




FIG. 17: Correlation function in the ground state at A = 1.2, 1.5. 
Straight lines are exponential laws e~ 1 ^. While for A = 1.5 the 
inverse correlation length = 0.0873 is very close to the value 
predicted by the Klein-Gordon model (A s /J = 0.0866), for A = 
1.2 it is considerably larger (£ _1 = 0.021, A s /J = 0.0048). 
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The constants C\ and C2 have been calculated in Ref [152]. 
Hence, in the whole planar phase (| A| < 1), the correlations 
exhibit critical behavior and fall off algebraically. 

A different situation has to be faced for A > 1. In 
the renormalization-group treatment backscattering terms be- 
come important. For A > 1, the sine-Gordon Hamiltonian, 



H. 



SG 



H 



LL 



2JA 
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dx cos(4</>(a:)) 



(B13) 



is the effective model. At the isotropic point, A = 1,K = 5, 
the cosine term is marginally relevant and leads to logarithmic 
corrections to the correlation function (B12). For Ising-like 
anisotropics, A > 1, < K < |, the cosine term is relevant 
- a spectral gap, A s , opens and the phase <j> becomes pinned 
at or ir/2. Hence, A = 1 marks a phase transition to an 
antiferromagnetically ordered phase with a finite asymptotic 
value of the spin-spin correlations, 



G zz (£) = (-1)V S 



(B14) 



The two degenerate ground states, corresponding to 4> = or 
7r/2, exhibit staggered magnetization, m s , of opposite signs, 



N ^ 



iy(S*) ~ (cos(20)) 



(B15) 



The spectral gap as well as the staggered magnetization are 
continuous in all derivatives in A - the phase transition is of 
Berezinskii-Kosterlitz-Thouless type [153, 154]. In Fig. 16, 
A s and m s are plotted as calculated from the Bethe ansatz 
[137, 138]. The energetically lowest excitations of the sine- 
Gordon model (B13) are solitons and antisolitons, which cre- 
ate kinks to antiferromagnetic domains with negativ (soli- 
tons), respectively positiv (antisolitons), sublattice magneti- 
zation m s . 

For sufficiently large anisotropies, where a semiclassical 
approximation becomes valid, the sine-Gordon model reduces 
essentially to the Klein-Gordon Hamiltonian, 
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As a result of the presence of the mass term in (B16), the 
connected correlation function, 

G Z7 \t) = ± Y,( S * S Ue) W) (S* +i ) , (B17) 



decays exponentially for large distances, 

G£*(*)~e-' /c , (B18) 
where the correlation length is given by the inverse gap, 

J 



A, 



(B19) 



In Fig. 17 the behavior (B18) is confirmed in the gapped state 
of the XXZ model by numerical simulations using imaginary- 
time evolution of the infinite-size matrix product state (see 



Ref. [155] or Section C for the description of this method). 
However, the relation (B19) is only valid for sufficiently large 
gaps. 

In order to avoid dealing with the complicated structure of 
the antiferromagnetic states in the XXZ model, we introduce 
the spin-density-wave (SDW) state, 

|#!DW= n k4+^ c !+ji°}- ( B2 °) 

-7r/2<fe<ir/2 

The coefficients of the wave function are related to the gap 
parameter A s by 
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(B21) 



The correlation function calculated with the state (B20) re- 
produces the exponential decay (B18) with the correlation 
length inversely proportional to the gap (B19) - the spin- 
density wave provides a valid phenomenological description 
of antiferromagnetic states. For special values of parameters 
- namely at the Luther- Emery point [151] - the spin-density 
wave (B20) coincides with the exact ground state of the sine- 
Gordon model (B13). Varying the gap parameter A s from 
zero to infinity, the spin-density-wave state (B20) links the 
ground state of the XXZ model at A = (B5) with the Neel 
state, 



i^>N« = iuT...m) 

the ground state in the limit A — > 00. 



(B22) 



APPENDIX C: MATRIX PRODUCT ALGORITHM FOR 
TIME-DEPENDENT PROBLEMS IN THE 
THERMODYNAMICS LIMIT 

The concept of matrix product states (MPS) [156-159] as 
a generalization of valence-bond states [143, 160, 161] has 
been developed parallel in time with the density matrix renor- 
malization group (DMRG) algorithm [162, 163]. DMRG 
established quickly as one of the most powerful numerical 
approaches for solving (quasi) one-dimensional correlated 
many-body problems at equilibrium. Although DMRG was 
originally introduced as a real-space renormalization group, it 
can be understood as a variational optimization procedure in 
the space of matrix product states [164]. This identification of 
DMRG with MPS is especially useful for the implementation 
of the ideas of DMRG in the thermodynamic limit [155, 165] 
and for time-dependent problems [72-74, 166]. In the fol- 
lowing we present a formulation of a DMRG-like algorithm, 
which is most suitable for both time-dependent and infinite- 
size calculations. The procedure is identical to the infinite-size 
time-evolving block decimation algorithm iTEBD [155], ex- 
cept that different matrices, introduced in the context of static 
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DMRG (Ref. [165]), are used in order to improve the stabil- 
ity of the algorithm. Since neither the density matrix nor the 
renormalization group idea appears explicitly in this formu- 
lation, we refer to the algorithm as the matrix product state 
algorithm MPS or iMPS, if the infinite-size limit shall be em- 
phasized. Error analysis will be given for a specific case of a 
non-equilibrium problem in the thermodynamic limit, where 
we find that the behavior of the error can be considered identi- 
cal to the case of the time-dependent DMRG for finite lattices 
[76]. 

1. Matrix product states 

In order to construct a MPS, we consider a one-dimensional 
lattice model where the Hilbert space can be separated into left 
and right subspaces L l and R l+1 - L 1 including i as rightmost 
site, i + 1 being the leftmost site of R %+1 . Generally, a wave 
function can be written as 

w=Ei*a 4 )( A *w*r l ). < ci > 

af) 

where ) + )) are orthonormal basis vectors of the 
space L l (R l+1 ). Aj is called the bond center matrix of bond 
i and constructs the density matrix of the left and right sub- 
systems, p L ' = AjAj and p R ' + = AiA\ respectively. If each 
site is described by a set of local basis vectors |sj) of dimen- 
sion cfj (si = 0, . . . , d{ — 1), a state of the subspace can be 
expanded in terms of the local basis and the remaining sub- 
space, 

|$f) = ^|< _1 }(A Si Wki}- (C2) 

The orthonormality of the basis imposes on Af the left or- 
thonormalization constraint, 

J2A?A(=6 SS ,. (C3) 

s 

Equivalently, the state of the right subspace can be expanded 
by means of right orthonormalized matrices, 

i$f) = Ei^( B rwi< +1 > 5 (C4) 

J2BfBf=5 ss ,. (C5) 

An iterative expansion of an arbitrary state on a lattice of 
size TV is possible, providing a matrix-product expression of 
the state, 

|^)=Tr £ A? A? ... A s N N \s lS2 ...s N ). (C6) 

S1S2...SN 

In the limit of N — > oo, in the presence of translational 
symmetry, a MPS can be constructed as a periodic array of 
matrices. Choosing a 2-site unit cell for concreteness, i.e. 



Ai + 2 = Ai, the set of matrices Af, A|, Bf, Ai, and 
A2 provides full information about the wave function. For in- 
stance, it is possible to construct the major object to be manip- 
ulated in a MPS algorithm, the two-site center matrix (index 
i denotes the site type, which can be either 1 or 2 for odd or 
even i), 

Af = A?AiB?' +1 , (C7) 

which can be used to decompose the wave function into the 
local bases of sites i and i + 1, 

m = E \*t%m a 'uw)\** i+ *)- cc8) 

ass' [3 

Left and right orthonormalized matrices are related via the 
single-site center matrix, 

AI = Ai_iB? = Af Ai , (C9) 

and can be formally mapped to each other, 

Bl = At\A\ and Af = Af A" 1 . (CIO) 

Using the single-site center matrix, the calculation of observ- 
ables for a MPS representation is straightforward. For a local 
operator Of s acting on site i, 

(Oi) = E °f (Tr A^ S 'A S ) . (CI 1) 

ss* 

Similarly, introducing an iterative procedure, correlation func- 
tions can be calculated [165]. 

2. Schmidt decomposition 

The preceding introduction of MPS is completely general 
and, if infinite-dimensional matrices are allowed, any state can 
be formally expressed in terms of a matrix product. A class of 
valence-bond states [143, 156-159, 161] is indeed naturally 
formulated in terms of MPS. Also, product states are trivially 
represented as MPS. Matrix product states are however espe- 
cially powerful in combination with an approximative numer- 
ical algorithm, providing the optimal reduced basis set for re- 
placing a large or possibly infinite Hilbert space. The Schmidt 
decomposition, as described in the following, is the procedure 
which allows to select the most relevant basis states. 

If only a finite number of states m is supposed to be re- 
tained (in order to keep the dimension of the Hilbert space 
manageable for the computer), it can be shown [163] that a 
state approximates best the targeted state \ip) in the form 
(CI), if it is defined as the Schmidt decomposition of rank m 
(site indices are omitted), 

m 

i^) = Ei |> «) a «i |> «)- ( C12 > 

a=l 

The Schmidt coefficients, X a , are the dominating eigenvalues 
of the singular value decomposition, 

A aP = E c / a 7 A 7 y; 7 , x\ > \\ > . . . , (ci3) 

7 
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satisfying J2 a = 1- The discarded weight, 

w=Y,*l> ( C14 ) 

corresponds to the mismatch, | \tp) — \tp) | = w, introduced by 
this truncation procedure. The new basis is given in terms of 
the Schmidt states, 

m) = j2 u ^) > 

a 

1*7) = E^l*»>- (C15) 

a 

In practice it is useful to set only an upper bound for m (rather 
than fixing a definite value) and instead define a threshold e 
such that only states for which > e are retained. The ap- 
plicability of this truncation procedure to a physical state de- 
pends on the characteristics of the Schmidt values or, equiv- 
alently, the spectrum of the density matrix. The more slowly 
the values A Q decay, the larger must be the number of retained 
states. A generic expression for the spectrum of the density 
matrix has been obtained for a critical theory [167], for prac- 
tical purposes it is however sufficient to consider the entangle- 
ment properties of the system to get the order of the necessary 
number of retained states. For instance, one can consider the 
entanglement entropy, 

S = Tr(p L ' R log 2 p L ' R ) = log 2 A* . (C16) 

a 

The fact that in one-dimensional equilibrium states the entan- 
glement entropy exhibits logarithmic dependence on the typi- 
cal length scale £ of the state [168] (£ corresponds to the corre- 
lation length or, at criticality, to the size of the system) guaran- 
tees an accurate description of a large class of wave functions 
using a finite number of states m oc £. Away from equilib- 
rium, however, the entanglement generally grows linearly in 
time [92] and a potentially exponential growth of m with time 
restricts the applicability of a MPS to short times. 

For a wave function represented at bond i by the two-site 
center matrix Af s , the Schmidt decomposition reads as fol- 
lows: Replacing in Eq. (CI) the contracted two-site center 
matrix with the bond center matrix, 

^diS+a,d iS '+p = (A* S ) a j3 , (C17) 

the Schmidt decomposition can be carried out as presented 
above. The matrices are updated retaining m Schmidt states 
(C15), 

{Al) a p — > U dl +a,f3 

(B? +1 ) a p -> V* !d2S2+ p (C18) 
(A,) a/ 3 — > 5 a p\ a , 

with Ai+i remaining unchanged. 

For the evaluation of correlation functions, additionally 
Aj + i or Bi are needed. Although the effect of loss of orthogo- 
nality is spurious when applying the direct inverse (CIO) as in 



the original iTEBD algorithm [155], especially in the case of 
real-time evolution, a procedure for recovering both left and 
right orthonormalized representations is needed for stabilizing 
the algorithm. The left (right) rotation of the matrices does the 
job. Starting from a single-site center matrix, Af = hi-\Bf, 
the left orthonormalized matrix and the rotated center matrix 
can be extracted from the singular value decomposition (C13) 
of the re-indexed matrix, A. a +ds,/3 = (A s ) a /3, 

{At) a0 = U a+ds . , (Af ) a/J = \ a Vp a . (C19) 

An iterative application of this procedure moves the center 
matrix through the lattice and brings all matrices into left or- 
thonormal form. An analogous left-moving iteration brings 
the matrices into the right orthonormalized form. In the pe- 
riodic iMPS a problem arises when the right-moving center 
matrix reaches the edge of the unit cell. Af does not in gen- 
eral coincide with the former Aj+i and repeating the itera- 
tions through the unit cell further changes the MPS. There ex- 
ist however schemes which solve this problem by introducing 
an additional transformation, after which the transfer opera- 
tor Af A"^ becomes equal to identity (see Refs. [165, 169] 
for detailed descriptions). 

3. Suzuki-Trotter decomposition 

In order to calculate the time evolution of a MPS, \ip(t)) = 
e~ lHt \ijjo), it is suitable to approximate the evolution opera- 
tor, e~ lHt , by a Suzuki-Trotter decomposition. This is pos- 
sible if the global operator H contains only nearest-neighbor 
bond terms H = ^ i/^i+i (e.g. the Heisenberg chain with 
Hi t i + i = SjSj+i). H can then be decomposed into even and 
odd parts H = Hi + H 2 , 

Hi = E • • • i H 2 = E H2j+i,2j+2 ■ (C20) 

i i 

The Suzuki-Trotter decomposition can be regarded as the first- 
order expansion of the evolution operator using the Baker- 
Hausdorff formula [170], 

e -iHt = ^-iHaig-iffitfjn + Q^n) , TlS = t . (C21) 

This approximation is improved in a second-order expansion, 

e -iHt = ( e - l H 1 S/2 e - l H 2 S e -iH 1 S/2- } n + ^3„J _ (C22) 

or, if higher accuracy is desired, using third- or higher-order 
expansions [171]. 

Since the components of the even (odd) part commute with 
each other, 

[H2j,2j+1, H2j',2j' + l] = [H 2 i+i :2 i+2, H2i' + 1.2i'+2] = 0, 

within the first-order Suzuki-Trotter decomposition the evo- 
lution operator can be broken down to a product of nearest- 
neighbor operators, 

e - iHt K I JJ e -^ i+1 ,2 i+ 2<5JJ e -iff2 3 ^ + l<5 I _ (C23 ) 



25 




Equivalent expressions hold for higher-order decompositions. 

The Suzuki-Trotter decomposition can also be used for cal- 
culating the ground state \ip) of a Hamiltonian using imagi- 
nary-time evolution, 



-tH 



(C24) 



where |^n) is some random initial state. In order to get reli- 
able results from this procedure, the Trotter slicing has to be 
reduced carefully during the imaginary-time evolution [155]. 



4. Update of an iMPS 

We consider now the application of a single factor of (C23) 
onto a MPS. For example, for the odd bond operator, U 

e~ 



, iH li2 S we naye 



(C25) 



After a subsequent singular value decomposition, retaining a 
finite number of states, the matrices can be updated, 



(C26) 



In the case of an iMPS of periodicity 2, the effect of the re- 
maining factors of e~ lHlS on the other odd bonds is identi- 
cal. Hence, the update (C26) corresponds to the action of the 



operator e 



-iHiS 



on the whole, infinitely extended wave func- 



tion. After the update one can recalculate Bf , A| by means 
of left- and right-moving iterations. More straightforward for 



preparing the application of the odd bond operator e 
the direct construction of the center matrix, 



-iH 2 S 



(C27) 



Since the inverse of A2 is required, a finite threshold e is nec- 
essary to guarantee the stability of this operation. 



5. Error analysis 

As an application of the iMPS method to a non-equilibrium 
problem, we study the quench problem in the XXZ chain, 



W)) 



-iHt 



\ijio), where |-0o) is the ground state of the 



XXZ Hamiltonian at a given value A = Ao, and H is char- 
acterized by an anisotropy parameter A. In this case the local 
basis consists of a spin-up and a spin-down states {| |), | f)}. 
This quench problem is analyzed in detail in section III. 

First we study the case where l^o) is the Neel state which 
has a trivial iMPS representation with m = 1, A\ = Bf = 
<5 s j, A2 — -Bf = 5 Since the z-projection of the total spin 
(St ot ) is conserved, the MPS can be resolved by this quantum 
number [165, 172]. The resulting speedup is about an order 
of magnitude in comparison with a simulation which exploits 
no symmetry. In the limit A = the numerical results can be 
checked against the exact solution (see section II). In Fig. 18 
the absolute deviation from the exact result, 



5m s {t) = |m,(t)- J (2J*)/2|, 



(C28) 



is plotted for different values of the number of retained states 
to, the threshold e and the Trotter slicing S. The evolution of 
the error can be clearly divided into two regimes by introduc- 



ing the runaway time t T 



For t < t r 



there is a 
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Jt 



FIG. 19: The absolute error (C29) in an iMPS simulation for S = 
0.005, e = 10 -17 for different values of m. Inset: the corresponding 
dynamics of m s (t). 

small error which does not depend on the value of to. In this 
case the error is dominated by the Trotter error, which grows 
at most linearly in time. However, for t > t runaway the error 
starts growing nearly exponentially. The approximately loga- 
rithmic dependence of t runaway on m (Fig. 18(a), inset) is in 
agreement with the linear growth of the entanglement entropy 
in the non-equilibrium problem - t runaway can be understood 
as the point where the chosen finite number of retained states 
is no more sufficient to represent the entanglement in the state. 
We note, however, that a strict relation between entanglement 
entropy and t runaway can not be rigorously established [76]. 

In order to reduce the Trotter error dominating at t < 
trunaway, one may choose smaller values of S. The thresh- 
old e has to be decreased as well. Otherwise, due to the 
increased number of updates, errors associated with the dis- 
carded weight at each step may accumulate. In Fig. 18(b) we 
plot two cases with e = 1CP 15 and 10~ 20 . In general, it is 
sufficient to reduce the threshold proportional to the Trotter 
slicing e oc St, 

If the Trotter slicing is chosen so that the resulting error is 



of the order of the accuracy goal, t runaway then sets the time 
window for the validity of the numerical results (in Fig. 18 
the accuracy goal in the absolute error is about 10 -6 ). We 
note that such behavior of the error in this time-dependent 
infinite-size MPS algorithm is identical to that of the finite- 
size DMRG algorithm [73]. 

If the exact solution is not known, t runaway can never- 
theless be determined by comparing curves from calculations 
with slightly different to. t runaway is the point where the dif- 
ference between them starts to grow significantly. Fig. 19 
illustrates this procedure with the results for a quench in the 
XXZ chain from A = 4 to A = 2, with S = 0.005 and 
e = 10 -17 (see also section III). Comparing the difference 
for various to, 

Sm s (t) = \m s (t) - m s (t) m=1400 | , (C29) 

where m s (£) m=1400 is the result for 1400 retained states, we 
find a behavior identical to the exactly solvable case of the XX 
chain - the curves for to < 1400 overlap completely with the 
one for to = 1400 up to t < t runaway and a difference can 
only be seen for t > t runaway . For to = 1400, t runaway is 
estimated in Fig. 19 by extrapolation of the values for to = 
600, 800, 1000. Again, the accuracy of the results for t < 
trunaway are dominated by the Suzuki-Trotter error which has 
to be estimated separately (here it is of the order of 10~ 7 ). 
In practice it is not mandatory to abort the calculation at the 
runaway time - from the rough behavior of 8m s (t) one can 
estimate that even for t < 10, the absolute error of the curve 
for m = 1400 is still of the order of 10~ 6 . 

The presented error analysis has been carried out for a lo- 
cal parameter in a specific setup. As long as non-equilibrium 
dynamics is concerned, this behavior is completely generic, 
although the runaway time and the Suzuki-Trotter error have 
to be determined for each case. Also, the error may depend 
on the observable under consideration - long-distance corre- 
lation functions may exhibit shorter runaway times than local 
observables. We would like to emphasize that error control, 
which imposes criteria on the wave function [173], is in gen- 
eral too strict, and the presented observable-based approach 
can considerably extend the accessible time window. 
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